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Numerical simulations of 15 orbits of an equal-mass binary black hole system are presented. 
Gravitational waveforms from these simulations, covering more than 30 cycles and ending about 
1.5 cycles before merger, are compared with those from quasi-circular zero-spin post-Newtonian 
(PN) formulae. The cumulative phase uncertainty of these comparisons is about 0.05 radians, 
dominated by effects arising from the small residual spins of the black holes and the small residual 
orbital eccentricity in the simulations. Matching numerical results to PN waveforms early in the 
run yields excellent agreement (within 0.05 radians) over the first ~ 15 cycles, thus validating the 
numerical simulation and establishing a regime where PN theory is accurate. In the last 15 cycles 
to merger, however, generic time-domain Taylor approximants build up phase differences of several 
radians. But, apparently by coincidence, one specific post-Newtonian approximant, TaylorT4 at 
3.5PN order, agrees much better with the numerical simulations, with accumulated phase differences 
of less than 0.05 radians over the 30-cycle waveform. Gravitational-wave amplitude comparisons 
are also done between numerical simulations and post-Newtonian, and the agreement depends on 
the post-Newtonian order of the amplitude expansion: the amplitude difference is about 6-7% for 
zeroth order and becomes smaller for increasing order. A newly derived 3.0PN amplitude correction 
improves agreement significantly (< 1% amplitude difference throughout most of the run, increasing 
to 4% near merger) over the previously known 2.5PN amplitude terms. 

PACS numbers: 04.25.D-, 04.25. dg, 04.25. Nx, 04.30.-w, 04.30. Db, 02.70.Hm 



I. INTRODUCTION 



The last two years have witnessed tremendous progress 
in simulations of black hole binaries, starting with 
the first stable simulation of orbiting and merging 
black holes [l|, development of the moving puncture 
method and rap id progress by other groups [1, 0, 

0, B H US Ellj [13, [U- Since then, an enormous amount 
of work has been done on the late inspiral and merger of 
black hole binaries, among them studies of the universal- 

J3i 



investigations into 



ity of the merger waveform s . , _ 

black hole kicks Tl|, [M [O, [iIdTzO, '21|, IH, [H, [H, [H, 
[H, [13, [Hi] and spin dynamjcs 29, 3Q, comparisons 
to post-Newtonian models 32, 33, 34| . and applications 
to gravitational wave data analysis |35l . [36l . |33| . 

Compared to the intense activity focusing on simu- 
lations close to merger, there have been relatively few 
simulations covering the inspiral phase. To date, only 
three simulations [H, [3^, [lo, [U, [4^ cover more than 
five orbits. Long inspiral simulations are challenging for 
a variety of reasons: First, the orbital period increases 
rapidly with separation, so that simulations must cover 
a significantly longer evolution time. In addition, the 
gravitational waveform must be extracted at larger ra- 
dius (and the simulation must therefore cover a larger 
spatial volume) because the gravitational wavelength is 
longer. Furthermore, gravitational wave data analysis 



requires small absolute accumulated phase uncertainties 
in the waveform, so the relative phase uncertainty of the 
simulation must be smaller. 

Gravitational wave detectors provide a major driv- 
ing force for numerical relativity (NR). The first gener- 
ation interferometric gravitational wave detectors, such 
as LIGO [4l,[43, GEO600 [H and VIRGO ^\^, are 
now operating at or near their design sensitivities. Fur- 
thermore, the advanced generation of detectors are en- 
tering their construction phases. This new generation 
of interferometers will improve detector sensitivity by a 
factor of ^ 10 and hence increase expected event rates 
by a factor of ~ 1000 4^. One of the most promising 
sources for these detectors is the inspiral and merger of 
binary black holes (BBHs) with masses mi ~ m2 ~ 10- 
20 Mq [i^. These systems are expected to have circu- 
larized long before their gravitational waves enter the 
sensitive frequency band of ground-based detectors [sOl ■ 

A detailed and accurate understanding of the gravita- 
tional waves radiated as the black holes spiral towards 
each other will be crucial not only to the initial de- 
tection of such sources, but also to maximize the in- 
formation that can be obtained from signals once they 
are observed. When the black holes are far apart, the 
gravitational waveform can be accurately computed us- 
ing a post-Newtonian (PN) expansion. As the holes ap- 
proach each other and their velocities increase, the post- 
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Newtonian expansion is expected to diverge from the true 
waveform. It is important to quantify any differences 
between theoretical waveforms and the true signals, as 
discrepancies will cause a reduction of search sensitiv- 
ity. Several techniques have been proposed to address 
the problem of the breakdown of the post-Newtonian ap- 
proximation [5TI , [5^ . [53| . but ultimately, the accuracy 
of the post-Newtonian waveforms used in binary black 
hole gravitational wave searches can only be established 
through comparisons with full numerical simulations. 

Unfortunately, comparing post-Newtonian approxima- 
tions to numerical simulations is not straightforward, the 
most obvious problem being the difficulty of producing 
long and sufficiently accurate numerical simulations as 
explained above. In addition, post-Newtonian waveforms 
typically assume circular orbits, and most astrophysi- 
cal binaries are expected to be on circular orbits late 
in their inspiral, so the orbital eccentricity within the 
numerical simulation must be sufficiently small^. An- 
other factor that complicates comparisons is the variety 
of post-Newtonian approximants available, from several 
straightforward Taylor expansions to more sophisticated 
Pade resummation techniques and the effective one-body 
approach (see e g. M, M,M,M,M,M,M,^, as 
well as Section Fill El below). While all post-Newtonian 
approximants of the same order should agree sufficiently 
early in the inspiral (when neglected higher-order terms 
are small) , they begin to disagree with each other during 
the late inspiral when the post-Newtonian approximation 
starts to break down — exactly the regime in which NR 
waveforms are becoming available. 

Finally, agreement (or disagreement) between NR and 
PN waveforms will also depend very sensitively on the 
precise protocol used to compare the waveforms. Are 
PN and NR waveforms matched early or late in the in- 
spiral? Is the matching done at a particular time, or is a 
least-squares fit performed over part (or all) of the wave- 
form? Does one compare frequencies uj{t) or phases (j){t)l 
Are comparisons presented as functions of time or of fre- 
quency? Up to which cutoff frequency does one compare 
PN with NR? 

Despite these difficulties, several comparisons between 
NR and PN have been done for the last few orbits of 
an equal-mass, non-spinning black hole binary. The 
first such study was done by Buonanno et al [32l | based 
on simulations performed by Pretorius ^ lasting some- 
what more than 4 orbits (~ 8 gravitational wave cycles). 
This comparison performs a least-squares fit over the full 
waveform, finds agreement between the numerical evo- 
lution and a particular post-Newtonian approximant (in 
our language TaylorTS 3.0/0.0^) and notes that another 



approximant (TaylorT4 3.5/0.0) will give similarly good 
agreement. However, as the authors note, this study is 
severely limited by numerical resolution, sizable initial 
eccentricity (~ 0.015), close initial separation of the black 
holes, and coordinate artifacts; for these reasons, the au- 
thors do not quantify the level of agreement. 

More recently. Baker et al. [H, [s^ performed simula- 
tions covering the last ~ 14 cycles before merger. These 
simulations have an orbital eccentricity ^ 0.008 [111 , forc- 
ing the authors to use a fitted smooth ( "de-eccentrized" ) 
gravitational wave phase to obtain a monotonically in- 
creasing gravitational wave frequency. Comparing to 
TaylorT4 3.5/2.5, they find agreement between numerical 
and post-Newtonian gravitational wave phase to within 
their numerical errors, which are about 2 radians. The 
authors also indicate that other post-Newtonian approxi- 
mants do not match their simulation as well as TaylorT4, 
but unfortunately, they do not mention whether any dis- 
agreement is significant (i.e., exceeding their numerical 
errors). Pan et. al [35| performed a more comprehensive 
analysis of the numerical waveforms computed by Preto- 
rius [32] and the Goddard group 38, 39] , confirming that 
TaylorT4 3.5/0.0 matches the numerical results best. 

The most accurate inspiral simulation to date was per- 
formed by the Jena group and presented in Husa et 
al. [42 1 and Hannam et al. [4l|. This simulation covers 



18 cycles before merger and has an orbital eccentricity of 
~ 0.0018 [6l|. Discarding the first two cycles which are 
contaminated by numerical noise, and terminating the 
comparison at a gravitational- wave frequency muj = 0.1 
(see Eq. ([T5|) for the precise definition) their comparison 
extends over 13 cycles. We discuss the results of Ref. [4l| 
in more detail in Sec. IVI A II 

This paper presents a new inspiral simulation of a non- 
spinning equal mass black hole binary. This new simu- 
lation more than doubles the evolution time of the sim- 
ulations in Refs. [H, [H, |4l|, resulting in a wave- 
form with 30 gravitational wave cycles, ending '--^ 1.5 cy- 
cles before merger, and improves numerical truncation 
errors by one to two orders of magnitude over those in 
Refs. [H, Hi], \^ . The orbital eccentricity of our sim- 
ulations is ~ 6 X 10^^; this low eccentricity is achieved 
using refinements of techniques described in [iO]. We 
present a detailed analysis of various effects which might 
influence our comparisons to post-Newtonian waveforms 
for non-spinning black hole binaries on circular orbits. 
These effects result in an uncertainty of ~ 0.05 radians 
out of the accumulated ^ 200 radians. Perhaps sur- 
prisingly, the largest uncertainty arises from the residual 
orbital eccentricity, despite its tiny value. The second 
largest effect arises due to a potential residual spin on 
the black holes, which we bound by \S\/M^^^ < 5 x 10"'*. 



^ Unfortunately, this circularization occurs on extremely long time 
scales |50( . thousands of orbits, making it impossible to run 
the numerical simulation long enough to radiate the eccentric- 
ity away. 

^ We identify post-Newtonian approximants with three pieces of 



information: the label introduced by [sSl for how the orbital 
phase is evolved; the PN order to which the orbital phase is 
computed; and the PN order that the amplitude of the waveform 
is computed. See Sec. IIII El for more details. 
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We compare the numerical waveforms with four 
different time-domain post-Newtonian Tayfor- 
approximants [s^, S, and we match PN and 

NR waveforms at a specific time during the inspiral. We 
explore the effects of varying this matching time. When 
matching ~ 9 cycles after the start of our evolution, all 
post-Newtonian approximants of 3.0PN and 3.5PN order 
in orbital phase agree with our simulation to within 
^ 0.03 radians over the first 15 cycles. This agreement 
is better than the combined uncertainties of the com- 
parison, thus validating our simulations in a regime 
where the 3.5PN truncation error of post-Newtonian 
theory is comparable to the accuracy of our simulations. 
Lower order post-Newtonian approximants (2.0PN and 
2.5PN order), however, accumulate a significant phase 
difference of 0.2 radians over this region. 

Extending the comparison toward merger (as well as 
when matching closer to merger), wc find, not surpris- 
ingly, that the agreement between PN and NR at late 
times depends stron gly on exactly what post-Newtonian 
approximant we use [54l . Issj . Typical accumulated phase 
differences are on the order of radians at frequency muj — 
0.1. One particular post-Newtonian approximant, Tay- 
lorT4 at 3.5PN order in phase, agrees with our NR wave- 
forms far better than the other approximants, the agree- 
ment being within the phase uncertainty of the compar- 
ison (0.05 radians) until after the gravitational wave fre- 
quency passes muj = 0.1 (about 3.5 cycles before merger). 
It remains to be seen whether this agreement is fun- 
damental or accidental, and whether it applies to more 
complicated situations (e.g. unequal masses, nontrivial 
spins). 

We also compare the post-Newtonian gravitational 
wave amplitude to the numerical amplitude, where we 
estimate the uncertainty of this comparison to be about 
0.5%. Restricted waveforms (i.e., OPN order in the am- 
plitude expansion) are found to disagree with the nu- 
merical amplitudes by 6-7%. An amplitude expansion of 
order 2PN shows significantly better agreement than the 
expansion at order 2.5PN. A newly derived 3PN ampli- 
tude (6^ is found to give much better agreement than 
the 2. OPN amplitude. 

This paper is organized as follows: Section HIl discusses 
our numerical techniques. In particular, we describe how 
we construct binary black hole initial data, evolve these 
data for 15 orbits, extract gravitational wave information 
from the evolution, and produce a gravitational waveform 
as seen by an observer at infinity. Section Hill details the 
generation of post-Newtonian waveforms, including de- 
tails of how we produce the four approximants that we 
compare against NR. We describe our procedure for com- 
paring NR and PN waveforms in Sec. llVi and present a 
detailed study of various sources of uncertainty in Sec.lVl 
The comparisons between NR and PN are presented in 
Section IVII This section is split into two parts: First, 
we compare each PN approximant separately with the 
numerical simulation. Subsequently, we show some addi- 
tional figures which facilitate cross-comparisons between 



the different PN approximants. Finally, we present some 
concluding remarks in Section rVIII The impatient reader 
primarily interested in NR-PN comparisons may wish to 
proceed directly to Table [TTTl summarizing the uncertain- 
ties of our comparisons, and then continue to Sec. IVII 
starting with Fig. [121 

II. GENERATION OF NUMERICAL 
WAVEFORMS 

In order to do a quantitative comparison between nu- 
merical and post-Newtonian waveforms, it is important 
to have a code capable of starting the black holes far 
enough apart to be in a regime where we strongly believe 
the post-Newtonian approximation is valid, track the or- 
bital phase extremely accurately, and do so efficiently so 
the simulation can be completed in a reasonable amount 
of time. Furthermore, the gravitational waves from such 
a simulation must be extracted in such a manner that 
preserves the accuracy of the simulation and predicts the 
waveform as seen by a distant observer, so a comparison 
with the post-Newtonian waveform can be made. In this 
section we describe the techniques we use to do this, as 
well as the results of a simulation starting more than 15 
orbits prior to merger. 

When discussing numerical solutions of Einstein's 
equations, we write all dimensioned quantities in terms 
of some mass scale m, which we choose to be the sum of 
the irreducible masses of the two black holes in the initial 
data: 

m = Afi„a + Mi„,2. (1) 

The irreducible mass of a single hole is defined as 

Mi„ = ^A/16n, (2) 

where A is the surface area of the event horizon; in prac- 
tice we take A to be the surface area of the apparent 
horizon. More generally, it is more appropriate to use 
the Christodoulou mass of each black hole, 

M|H = M2,-f-^, (3) 

in- 

instead of the irreducible mass. Here S is the spin of the 
hole. However, for the case considered in this paper, the 
spins are sufficiently small that there is little difference 
between Mbh and Mi„. 

A. Initial data 

Initial data are constructed within the conformal thin 
sandwich formalism [H, [63| using a pseudo-spectral el- 
liptic solver [65|. We employ quasi-equilibrium bound- 
ary conditions|66l. [67| on spherical excision boundaries, 
choose conformal flatness and maximal slicing, and use 
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FIG. 1: Proper separation (top panel) and its time derivative 
(lower panel) versus time for short evolutions of the d = 30 ini- 
tial data sets 30a, 30b, and 30c (see Table |T]). These three data 
sets represent zero through two iterations of our eccentricity- 
reduction procedure. The orbital eccentricity is reduced sig- 
nificantly by each iteration. 



Eq. (33a) of Ref. '68'] as the lapse boundary condition. 
The spins of the black holes are made very small via an 
appropriate choice of the tangential shift at the excision 
surfaces, as described in (GSj . 

As the most accurate post-Neviftonian waveforms avail- 
able assume adiabatic inspiral of quasi-circular orbits, it 
is desirable to reduce the eccentricity of the numerical 
data as much as possible. Using techniques developed 
in [i^l, each black hole is allowed to have a nonzero ini- 
tial velocity component towards the other hole. This 
small velocity component Vr and the initial orbital angu- 
lar velocity Qq are then fine-tuned in order to produce 
an orbit with very small orbital eccentricity^. We have 
improved our eccentricity-reduction procedure since the 
version described in [40| . so we summarize our new iter- 
ative procedure here: 

We start with a quasi-circular (i.e., Vr = 0) initial data 
set at coordinate separation d = 30, where is de- 
termined by equating Komar mass with Arnowitt-Deser- 
Misner (ADM) mass ^62,]. We then evolve these data for 
about 1.5 orbits, corresponding to a time t/m fa 600. 
From this short evolution, we measure the proper sepa- 
ration s between the horizons by integration along the 



^ An alternative method of producing low-eccentricity initial data, 
based on post-Newtonian ideas, is developed in [611 . While that 
technique is computationally more efficient than ours, it merely 
reduces orbital eccentricity by a factor of ~ 5 relative to quasi- 
circular initial data, which is insufficient for the comparisons pre- 
sented here, (cf Sec. lVE2l l. 



coordinate axis connecting the centers of the black holes. 
We fit the time derivative ds/dt in the interval 100 < 
t/m ^ 600 to the function 



ds 
lit 



Aq + Alt + Bcos{ujt + if), 



(4) 



where we vary all five parameters Aq, Ai,B,uj and Lp to 
achieve the best fit. The desired smooth inspiral is rep- 
resented by the part Aq Ait\ the term B cos{u)t + ip) 
corresponds to oscillations caused by orbital eccentricity. 

For a Newtonian orbit with radial velocity B cos{ujt + 
Lp) at initial separation sq, it is straightforward to deter- 
mine the changes to the orbital frequency and the radial 
velocity which make the orbit perfectly circular, namely 



fJo VLq 



B sin 

2so ' 
BcosLp 



(5) 
(6) 



For Newtonian gravity, Eq. ^ will of course result in 
a circular orbit with Vr = 0. In General Relativity, 
and Vr will be different from their Newtonian values, for 
instance < to account for the inspiral of the two 
black holes. Nevertheless, we assume that small pertur- 
bations around the zero-eccentricity inspiral trajectory 
behave similarly to small perturbations around a Newto- 
nian circular orbit. Therefore, we apply the same formu- 
lae, Eqs. (O and to obtain improved values for Slg 
and Vr for the black hole binary, where sq is the initial 
proper separation between the horizons. We then use the 
new values of Qq and Vr to construct a new initial data 
set, again evolve for two orbits, fit to Eq. ([4]), and update 
flo and Vr- We continue iterating this procedure until the 
eccentricity is sufficiently small. 

We estimate the eccentricity for each iteration from 
the fit to Eq. g]) using the formula 



Sds/dt 



B 



(7) 



which is valid in Newtonian gravity for small eccentric- 
ities. Successive iterations of this procedure are illus- 
trated in Fig. [T] and yield the initial data sets 30a, 30b, 
and 30c summarized in TablelH Eccentricity decreases by 
roughly a factor of 10 in each iteration, with 30c having 
^ds/dt ~ 5 X 10^^. The evolutions used during eccentric- 
ity reduction need not be very accurate and need to run 
only for a short time, t ^ 600to. One iteration of this pro- 
cedure at our second lowest resolution requires about 250 
CPU-hours. For completeness, Table H] also lists parame- 
ters for initial data at smaller separation; these data will 
be used for consistency checks below. Apart from these 
consistency checks, the remainder of this paper will focus 
exclusively on evolutions of the low-eccentricity initial 
data set 30c. 
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TABLE I: Summary of the initial data sets used in this paper. The first block of numbers (d, Qo, fr, and Vr) represent raw 
parameters entering the construction of the initial data. The second block gives some properties of each initial data set: m 
denotes the sum of the irreducible masses, Madm and Jadm the ADM energy and angular momentum, and sq the initial proper 
separation between the horizons. The last column lists the eccentricity computed from Eq. ((7}. The initial data set 30c is used 
for all evolutions (except for consistency checks) described in this paper. 



Name 


d 




fr 


Vr X 10" 




MAUM/m 


JAum/m? 


So/m 


^ds/dt 


30a 


30 


0.0080108 


0.939561 


0.00 


0.01664793 


0.992333 


1.0857 


17.37 


1.0 X 10 


-2 


30b 


30 


0.0080389 


0.939561 


-4.90 


0.0167054 


0.992400 


1.0897 


17.37 


6.5 X 10 


-4 


30c 


30 


0.0080401 


0.939561 


-4.26 


0.0167081 


0.992402 


1.0898 


17.37 


5 X 10 


5 


24a 


24 


0.0110496 


0.92373 


-8.29 


0.0231947 


0.990759 


1.0045 


14.15 


1.1 X 10 


-3 


24b 


24 


0.0110506 


0.923739 


-8.44 


0.0231967 


0.990767 


1.0049 


14.15 


1.5 X 10 


-4 



B. Evolution of the inspiral phase 

The Einstein evolution equations are solved with the 
pseudo-spectral evolution code described in Ref. 0] ■ This 
code evolves a first-order representation [6^ of the gen- 
eralized harmonic system [tO, ItT], We handle the 
singularities by excising the black hole interiors from our 
grid. Our outer boundary conditions [H, [tI, [tI] are de- 
signed to p revent the influx of unphysical constraint vio- 
lations [7g.l76ll77l. l7i . l79ll8Ol . l8ll and undesired incoming 
gravitational radiation |82| , while allowing the outgoing 
gravitational radiation to pass freely through the bound- 
ary. 

The code uses a fairly complicated domain decompo- 
sition to achieve maximum efficiency. Each black hole is 
surrounded by several (typically six) concentric spher- 
ical shells, with the inner boundary of the innermost 
shell (the excision boundary) just inside the horizon. A 
structure of touching cylinders (typically 34 of them) sur- 
rounds these shells, with axes along the line between the 
two black holes. The outermost shell around each black 
hole overlaps the cylinders. The outermost cylinders 
overlap a set of outer spherical shells, centered at the ori- 
gin, which extend to large outer radius. External bound- 
ary conditions are imposed only on the outer surface of 
the largest outer spherical shell. We vary the location of 
the outer boundary by adding more shells at the outer 
edge. Since all outer shells have the same angular resolu- 
tion, the cost of placing the outer boundary farther away 
(at full resolution) increases only linearly with the ra- 
dius of the boundary. External boundary conditions are 
enforced using the method of Bjorhus [83], while inter- 
domain boundary conditions are enforced with a penalty 
method 

We employ the dual-frame method described in 
Ref. [7!]: we solve the equations in an 'inertial frame' 
that is asymptotically Minkowski, but our domain de- 
composition is fixed in a 'comoving frame' that rotates 
with respect to the inertial frame and also shrinks with 
respect to the inertial frame as the holes approach each 
other. The positions of the holes are fixed in the comov- 
ing frame; we account for the motion of the holes by dy- 
namically adjusting the coordinate mapping between the 
two frames. Note that the comoving frame is referenced 



only internally in the code as a means of treating mov- 
ing holes with a fixed domain. Therefore all coordinate 
quantities (e.g. black hole trajectories, wave-extraction 
radii) mentioned in this paper are inertial-frame values 
unless explicitly stated otherwise. 

One side effect of our dual frame system is that the 
outer boundary of our domain (which is fixed in the co- 
moving frame) moves inward with time as observed in 
the inertial frame. This is because the comoving frame 
shrinks with respect to the inertial frame to follow the 
motion of the holes. In Refs. 0, S^l the inertial frame 
coordinate radius r (with respect to the center of mass) 
and the comoving coordinate radius r' are related by a 
simple scaling 



r — a{t)r' . 



(8) 



The expansion parameter a{t) is initially set to unity 
and decreases dynamically as the holes approach each 
other, so that the comoving-frame coordinate distance 
between the holes remains constant. The outer boundary 
of the computational grid is at a fixed comoving radius 
^bdry which is mapped to the inertial coordinate radius 
-Rbdry(^) = a{t)R'^^^^^. Bccausc we wish to accurately 
compute the gravitational radiation as measured far from 
the holes, it is desirable to have a moderately large outer 
boundary (i?bdry(t) ^ 200rn) throughout the run. For 
the linear mapping, Eq. this requires a very dis- 
tant outer boundary early in the run, i?bdry(0) — 1000m. 
Computationally this is not very expensive. However, 
the initial junk radiation contaminates the evolutions for 
a time interval proportional to the light-crossing time to 
the outer boundary, and for i?bdry(0) — 1000m it would 
be necessary to discard a significant portion of the evo- 
lution. 

We therefore use the mapping 



i{t) + (1 - a{t)) -2 



(9) 



for some constant R'q which is chosen to be roughly the 
radius of the outer boundary in comoving coordinates. 
This mapping has the following properties: (1) At the 
initial time t — Q, the map reduces to the identity map 
because a(0) = 1. Thus we do not need to re-map our 



TABLE 11: Overview of low-eccentricity simulations discussed 
in this paper. -Rbdry is the initial coordinate radius of the 
outer boundary; this radius changes during the evolution ac- 
cording to the choice of "radial map" between inertial and 
comoving coordinates. The last column lists the different res- 
olutions run for each evolution, N6 being highest resolution. 
Evolution 30c-l/N6 forms the basis of our post-Newtonian 
comparisons, and is used in all figures unless noted otherwise. 



Name 


ID 


^''orbits 


-Rbdry radial map 


resolutions 


30c-l 


30c 


15.6 


462m Eq. ^ 


Nl, N2, N6 


30c-2 


30c 


15.6 


722m Eq. © 


N2, N4, N6 


30c-3 


30c 


15.6 


202m Eq. JH} 


N2, N3, . . . , N6 


24b- 1 


24b 


8.3 


160m Eq. JS} 


N2, N3, N4 




FIG. 2: Spacetime diagram showing the spacetime volume 
simulated by the numerical evolutions listed in Tab. [Ill The 
magnified view in the right panel shows how the gravitational 
waves are escorted to our extraction radii (see Sec. Ill C|) after 
the simulation in the center has already crashed at t ~ 3930m, 
and after the estimated time of the black hole merger, which 
is indicated by the circle. The thin diagonal lines are lines 
of constant t — r*; each corresponds to a retarded time at 
which the gravitational wave frequency ui at infinity assumes 
a particular value. 



initial data before evolving. (2) For small radii (i.e., at 
the locations of the black holes), the map reduces to the 
linear map, r = a{t)r' + 0{r'^). This allows use of the 
control system without modifications. (3) The moving 
radius r' = R'q is mapped to a constant inertial radius: 
r(i?Q) = Rq. This allows us to keep the inertial radius 
of the outer boundary constant (or nearly constant^) in 
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-4 
x/m 



FIG. 3: Coordinate trajectories of the centers of the black 
holes. The small circles/ellipsoids show the apparent horizons 
at the initial time and at the time when the simulation ends 
and wave escorting begins. The inset shows an enlargement 
of the dashed box. 



time rather than shrinking rapidly. 

In total, we have run three evolutions of the 30c ini- 
tial data set; these use different combinations of outer 
boundary radius and radial mapping between inertial and 
moving coordinates. Some properties of these evolutions 
are summarized in Table [III We also performed exten- 
sive convergence testing, running the same evolution on 
up to six distinct resolutions, Nl to N6. The coarsest 
resolution 30c-l/Nl uses approximately 41'^ grid points 
(summing all grid points in all the subdomains) , while the 
most accurate evolution, 30c-l/N6, uses about 67'^ grid 
points. The run 30c-l/N2 required about 2,500 CPU- 
hours and run 30c-l/N6 about 19,000, where our simula- 
tions do not take advantage of symmetries. The distance 
to the outer boundary is adjusted by adding or removing 
outer spherical shells to an otherwise unmodified domain- 
decomposition. Run 30c-l has 20 such outer spherical 
shells, while 30c-2 utilizes 32 and 30c-3 only 8. Thus, the 
total number of grid points varies slightly between runs, 
e.g. about 71'^ for 30c-2/N6. Figure [2] indicates the dif- 
ferent behavior of the outer boundary location for these 
three evolutions. 

For all of the evolutions 30c-l/2/3, the coordinate tra- 
jectories of the centers of the apparent horizons appear 
as in Fig. [3] The regular inspiral pattern without notice- 



* In practice, we choose R'q somewhat larger than the outer bound- 
ary, so that the outer boundary of the computational domain 



slowly contracts in inertial coordinates. This makes the zero- 
speed characteristic fields outgoing there, avoiding the need to 
impose boundary conditions on those fields. 
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able oscillations once again indicates that our evolutions 
indeed have very low eccentricity. 

Figure S] demonstrates the convergence of the black 
hole mass m(t) with spatial resolution for run 30c-l. The 
mass m(t) is computed as the sum of the irreducible 
masses of both black holes, as defined in Eq. ([2]). At 
the highest resolution, m{t) deviates by only a few parts 
in 10^ from its initial value m. 

Our apparent horizon finder works by expanding the 
radius of the apparent horizon as a series in spherical 
harmonics up to some order L. We utilize the fast flow 
methods developed by Gundlach [s^ to determine the 
expansion coefficients; these are significantly faster than 
our earlier minimization algorithms [stI . [s^ . The appar- 
ent horizon is almost spherical during the inspiral, so that 
the expansion in L converges exceedingly fast: L = 8 re- 
sults in a relative error of the irreducible mass of better 
than 10~*. The distortion of the horizons becomes more 
pronounced toward the end of the evolution when the 
black holes approach each other rapidly. This results in 
an error of 10~^ in the L — 8 apparent horizon expansion 
for the last 10m of the evolution. 

We also measure the quasi-local spin using coordinate 
rotation vect ors p rojected into the apparent horizon sur- 
faces [H, HO, mf. Only the z-component of the spin is 
non-zero (i.e., the spins are aligned with the orbital angu- 
lar momentum). The spin starts at S./M^.^ « -6 x 10^^ 
and increases slowly to —5 x 10""' during the evolution, 
where the minus sign indicates that the black hole spin is 
anti-aligned with the orbital angular momentum. Thus it 
appears the black hole's spins move further away from the 
corotational state. We believe this effect is caused by the 
use of coordinate rotation vectors when calculating the 
quasi-local spin, rather than more sophisticated approxi- 
mate Killing vectors [qIjIoI,!!^. Preliminary results with 
approximate Killing vectors find the initial spin to be less 
than 10~^, and slowly increasing during the evolution to 
a final value of 2 x 10~^ at the end of the comparison in- 
terval to post-Newtonian theory. Given the preliminary 
character of these results, we will take here the conserva- 
tive bound |S|/Af^j, < 5 x 10""* obtained from coordinate 
rotation vectors. 



C. Escorting gravitational waves 

The simulation presented in Fig.[3]stops when the hori- 
zons of the black holes become too distorted just before 
merger. At this time, most of the domain (all regions 
except for the immediate vicinity of the two holes) is 
still well resolved, and the spacetime contains gravita- 
tional radiation that has not yet propagated out to the 
large radii where we perform wave extraction. So in- 
stead of losing this information, which consists of sev- 
eral gravitational-wave cycles, we evolve only the outer 
portions of our grid beyond the time at which the code 
crashes in the center, effectively 'escorting' the radiation 
out to the extraction radii. 
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FIG. 4: Deviation of total irreducible mass m(f) = 2Mirr(t) 
from its value in the initial data. Plotted are the six different 
resolutions of run 30a- 1. 

To do this, we first stop the evolution shortly be- 
fore it crashes, and we introduce a new spherical exci- 
sion boundary that surrounds both black holes and has 
a radius of roughly three times the black hole separa- 
tion. This new excision boundary moves radially out- 
ward at slightly faster than the speed of light so that it 
is causally disconnected from the interior region where 
the code is crashing, and so that no boundary conditions 
are required on this boundary. We then continue the 
evolution on the truncated spherical-shell domain that 
extends from the new excision boundary to the outer 
boundary. To move both boundaries appropriately, we 
employ a new radial coordinate mapping 

r = ^(t)r(r') +B(t), (10) 

where r(r') is given by Eq. ([5]). The functions A(t) 
and B(t) are chosen to satisfy three criteria: First, the 
inner boundary of the spherical shell moves outward 
with coordinate speed of unity, which turns out to be 
slightly superluminal. Second, the outer boundary loca- 
tion -Rbdry(0 continuous first and second time deriva- 
tives at the time we transition to the truncated domain. 
And finally, the outer boundary location i?bdry(i) ap- 
proaches some fixed value at late times. The right panel 
of Fig. [2] shows the motion of the inner and outer radii 
for evolutions 30c- 1 and 30c-2 (we did not perform wave 
escorting for 30c-3). For 30c-l, wave escorting extends 
the evolution for an additional time 220m beyond the 
point at which the simulation stops in the center. 

Figure [H] shows the gravitational waveform extracted 
at inertial coordinate radius R ~ 2A0m for the run 30c- 
1. The brown vertical line indicates the time when wave 
escorting starts. Wave escorting allows us to extract an- 
other 4 cycles of gravitational waves. When comput- 
ing the gravitational wave strain h(t) from the Newman- 
Penrose scalar ^4 (see Eq. pT|) below), one must choose 
integration constants during the time integration. These 
integration constants were chosen such that h{t) has zero 
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FIG. 5: Gravitational waveform extracted at r = 240m. From 
top panel to bottom: The real part of the (2, 2) component 
of r^i; the gravitational wave strain, obtained by two time 
integrals of Re(r>I'4); the frequency of the gravitational wave, 
Eq. (I15|l : the gravitational wavelength, A = 2n/ij. The verti- 
cal brown line a.t t ~ 3930m indicates the time when "wave 
escorting" starts. 




t/m 



FIG. 6: Constraint violations of run 30c- 1. The top panel 
shows the norm of all constraints, normalized by the 
norm of the spatial gradients of all dynamical fields. The 
bottom panel shows the same data, but without the normal- 
ization factor. Norms are taken only in the regions outside 
apparent horizons. 



average and first moment [40| , which is is sufficiently ac- 
curate for the illustrative Fig. El To avoid errors caused 
by the choice of integration constants, the comparison 
to post-Newtonian waveforms below is based entirely on 

\E'4. 

In the lower two panels of Fig. [5] there is a signifi- 
cant amount of noise near the beginning of the run, at 
t < 250to. This noise is barely evident in the top panel 
of Fig. [5] as well. The noise is a manifestation of 'junk ra- 
diation', a pulse of radiation often seen at the beginning 
of numerical relativity simulations, and is caused by the 
initial data not being precisely a snapshot of an evolution 
that has been running for a long time. Among the effects 
that produce junk radiation are incorrect initial distor- 
tions of the individual holes, so that each hole radiates 
as it relaxes to its correct quasi-equilibrium shape. 

Our evolution code docs not explicitly enforce either 
the Einstein constraints or the secondary constraints that 
arise from writing the system in first-order form. There- 
fore, examining how well these constraints are satisfied 
provides a useful consistency check. Figure [S] shows the 
constraint violations for run 30c- 1. The top panel shows 
the norm of all the constraint fields of our first or- 
der generalized harmonic system, normalized by the 
norm of the spatial gradients of the dynamical fields (see 
Eq. (71) of Ref. ^). The bottom panel shows the same 
quantity, but without the normalization factor {i.e., just 
the numerator of Eq. (71) of Ref. [1^). The norms 
are taken over the entire computational volume that lies 
outside of apparent horizons. At early times, t < 500m, 



the constraints converge rather slowly with resolution be- 
cause the junk radiation contains high frequencies. Con- 
vergence is more rapid during the smooth inspiral phase, 
after the junk radiation has exited through the outer 
boundary. The constraints increase around t ^ 3900m 
as the code begins to fail near the two merging holes, 
but then the constraints decrease again after the failing 
region is excised for wave escorting. The normalized con- 
straint violations are less than 10~^ until just before the 
peak (which occurs at t = 3930m for all but the lowest 
resolutions). The size of the peak causes some concern 
that the waveforms at late times may be contaminated 
by constraint violations to a non-negligible degree. How- 
ever, near the peak, the constraint violations are large 
only in the inner regions of the domain near the black 
holes (note that the curves in Fig. [5] decrease by two or- 
ders of magnitude immediately after these inner regions 
are excised at t = 3930m). Because all constraint quan- 
tities propagate at the speed of light or slower for the 
formulation of Einstein's equations that we use, any in- 
fiucnce that the constraint peak has on the extracted 
waveform occurs after the constraint violations have had 
time to propagate out to the wave extraction zone. This 
is very late in the waveform, well after the gravitational 
wave frequency reaches mu> = 0.1, as can be seen from 
the right panel of the spacetime diagram in Fig. [21 
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D. Waveform extraction 

Gravitational waves are extracted using the Newman- 
Penrose scalar ^"4, using the same procedure as in [40| . 
To summarize, given a spatial hypersurface with timelike 
unit normal n'^, and given a spatial unit vector r'^ in the 
direction of wave propagation, the standard definition 
of \l/4 is the following component of the Weyl curvature 
tensor, 

*4 = -Ca^/5,rrm"m'3, (11) 

where i?^ = -^{'^^ a-nd is a complex null vector 

(satisfying to^to^ = 1) that is orthogonal to and . 
Here an overbar denotes complex conjugation. 

For (perturbations of) flat spacetime, 'I'4 is typically 
evaluated on coordinate spheres, and in this case the 
usual choices for ri^, r'' and are 

(12a) 

where (r, 9, (ji) denote the standard spherical coordinates. 
With this choice, ^4 can be expanded in terms of spin- 
weighted spherical harmonics of weight —2: 

*4(i, r, e,cf>)^Y. *4"(i, r) -2Yi^{0, cj)), (13) 

Im 

where the ^'4™ are expansion coefficients defined by this 
equation. 

For curved spacetime, there is considerable freedom in 
the choice of the vectors r^^ and m^, and different re- 
searchers have made different choices [1, IH, [9^ [9^ 113, 
[98. . ,99J that are all equivalent in the r ^ 00 limit. We 
choose these vectors by first picking an extraction two- 
surface £ that is a coordinate sphere (r^ = x'^+y'^ + z'^ us- 
ing the global asymptotically Cartesian coordinates em- 
ployed in our code) centered on the center of mass of the 
binary system, i.e. the point of symmetry. We choose r'' 
to be the outward-pointing spatial unit normal to £ (that 
is, we choose proportional to V^r and raise the index 
with the spatial metric). Then we choose according 
to Eq. (|12cp . using the standard spherical coordinates 9 
and (f> defined on these coordinate spheres. Finally we 
use Eqs. ini) and ([13]) to define the coefficients. 

Note that the m'^ vector used here is not exactly null 
nor exactly of unit magnitude at finite r. The resulting 
^4" at finite r will disagree with the waveforms observed 
at infinity. Our definition does, however, agree with the 
standard definition given in Eqs. pT|) - (fT3|) as r ^ 00. 
Because we extrapolate the extracted waves to find the 
asymptotic radiation field (see Section llIF[) . these effects 
should not play a role in our FN comparisons: Relative 



errors in ^'4 introduced by using the simple coordinate 
tetrad fall off like 1 /r, and thus should vanish after ex- 
trapolating to obtain the asymptotic behavior. While 
more careful treatme nt of the extraction method — such 
as those discussed in llOll Il02| | — may improve the 
quality of extrapolation and would be interesting to ex- 
plore in the future, the naive choice made here should be 
sufficient to ensure that the waveform after extrapolation 
is correct to the accuracy needed for these simulations. 

In this paper, we focus on the {l,m) = (2,2) mode. 
Following common practice (see e.g. [1,113), split the 
extracted waveform into real phase cj) and real amplitude 
A, defined by 

^f{r,t)=A{r,t)e-"^^''''\ (14) 
The gravitational-wave frequency is given by 

^ = i 

The minus sign in the definition of (f) is chosen so that the 
phase increases in time and m is positive. Equation (jl4p 
defines 4> only up to multiples of 27r. These multiples of 
27r are chosen to make (j) continuous through each evolu- 
tion, still leaving an overall multiple of 27r undetermined. 
We will consider only phase differences in this paper, so 
the choice of this overall phase offset is irrelevant. 

E. Convergence of extracted waveforms 

In this section we examine the convergence of the grav- 
itational waveforms extracted at fixed radius, without 
extrapolation to infinity. This allows us to study the 
behavior of our code without the complications of ex- 
trapolation. The extrapolation process and the resulting 
extrapolated waveforms are discussed in Sec. IIIFI 

The top panel of Fig. [7] shows the convergence of the 
gravitational wave phase (j) with numerical resolution for 
the run 30c-l. For this plot, the waveform is extracted at 
a fixed radius R = 77m. Each line shows the difference 
between </> computed at some particular resolution and 
(j) computed from our highest-resolution run 30c-l/N6. 
When subtracting results at different resolutions, no time 
or phase adjustment has been performed. The difference 
in between the two highest-resolution runs is smaller 
than 0.03 radians throughout the run, and it is smaller 
than 0.02 radians between t = IOOOto and the point at 
which muj = 0.1. 

At times before 1000m, the phase convergence of our 
simulation is limited to about 0.05 radians because of 
effects of junk radiation (described at the end of Sec- 
tion [irC])- The sharp pulse of junk radiation has com- 
paratively large numerical truncation error, and excites 
all characteristic modes at truncation-error level, includ- 
ing waves that propagate back toward the origin. Gen- 
eration of these secondary waves stops when the pulse of 
junk radiation leaves through the outer boundary (i.e., 
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FIG. 7: Convergence of the gravitational wave phase ex- 
tracted at radius R = 77m. All lines show differences with re- 
spect to our highest resolution run, 30c-l/N6. The top panel 
shows different resolutions of the same run 30c-l; no time or 
phase shifts have been performed. The bottom panel com- 
pares different runs, aligning the runs at mui = 0.1 by a time 
and phase shift. The thin vertical line indicates the time at 
which muj = 0.1 for 30c-l/N6. 



after one light-crossing time). Because we use the im- 
proved outer boundary conditions of Rinne et al. [t^ . 
there are no significant reflections when the junk radia- 
tion passes through the outer boundary. However, the 
waves produced before the junk radiation leaves remain 
in the computational domain for two additional light- 
crossing times, until they eventually leave through the 
outer boundary. 

The bottom panel of Fig. [7] shows phase comparisons 
between different waveforms after we perform a time shift 
and phase shift so that the waveforms agree at muj = 0.1. 
Our procedure for time shifting and phase shifting is the 
same as the shifting procedure we use to compare NR 
with PN waveforms (see Sec. lIVBp . so that the error 
estimates we extract from the bottom panel of Fig. [7] are 
relevant for our later NR-PN comparison. 

There arc three different types of comparisons shown 
in the bottom panel of Fig. [7l Phase differences between 
runs with the same initial data but with different outer 
boundary locations, phase differences between runs with 
different initial data, and phase differences between dif- 
ferent numerical resolutions of the same run (this last 
comparison is the same as what is shown in the top panel, 
except in the bottom panel the waveforms are time and 
phase shifted) . We will discuss all three of these in turn. 

First, we compare the phase difference of 30c-l/N6 
with runs that have different outer boundary locations. 
Run 30c-2 (with more distant outer boundary) agrees to 
within 0.002 radians with run 30c-l, but run 30c-3 (with 



closer outer boundary), has a much larger phase differ- 
ence with 30c- 1. We believe that this is because run 30c- 
3 has a very small ratio of outer boundary location to 
gravitational wavelength: R/X is about 1.1 for the first 
two-thirds of the run, and remains less than 2 for the 
entire run. 

We can explain the order of magnitude of these phase 
differences using the analysis of Buchman & Sarbach . 
Our outer boundary conditions are not perfectly absorb- 
ing, but instead they reflect some fraction of the outgoing 
radiation.^ The ratio of the amplitude of curvature per- 
turbations (i.e. ^'4) of the reflected wave to that of the 
outgoing wave is 



2(27r)4 



(16) 



The incoming reflected waves grow like 1 /r as they travel 
inward just like the outgoing waves decrease by l/r as 
they propagate outward. Therefore, the ratio of ampli- 
tudes of incoming and outgoing waves will have approx- 
imately the same value, q, at smaller radii, and we as- 
sume for the sake of this rough argument that this ratio 
remains equal to q even in the vicinity of the black holes 
(where it is no longer technically meaningful to talk about 
'radiation'). Now consider the second time derivative of 
the gravitational wave phase, (/>; this is nonzero only be- 
cause of gravitational wave emission, so (j) is proportional 
to some power of the outgoing wave amplitude. To get 
the correct power, we can use Eq. (|45p to find x x^, 
so Eq. yields (I) ~ x^^/^ (we assume gravitational 
wave phase is twice the orbital phase). The amplitude 
of 5*4 scales like x^, so (j) ^ A^^^^. Let us assume for 
the sake of this rough error estimate that the change in 
(f> due to the ingoing reflected wave scales similarly with 
amplitude, ~ A^^^^, where A = qA is the amplitude 
of the reflected ingoing wave. Therefore the unphysical 
gravitational-wave force acting back on the system due 
to boundary reflections will cause fractional errors in the 
second derivative of the phase of about q^^^^. That is, 
the error Scj) caused by the improper boundary condition 
will be given by 



,11/ 



(17) 



Integrating this yields 6<j) — q '^(j), where (f> is the total 
gravitational wave phase accumulated during the evolu- 
tion. For 30c-3, \/R - 0.9, so g ~ 6 x 10""*, which yields 
6(1) ^ 0.08 radians for an accumulated gravitational wave 
phase of about 200 radians. This rough estimate agrees 
in order of magnitude with the phase difference between 



^ However, in a comparison of various boundary conditions [74l |. 
the boundary conditions we use produced smaller reflections than 
other boundary conditions commonly used in numerical relativ- 
ity 
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30c-3 and 30c-l as shown in the bottom panel of Fig- 
ure [71 The run 30c- 1 has an outer boundary about 2.5 
farther away, reducing the reflection coefficient by a fac- 
tor 2.5^ « 40, so for 30c-l this estimate of the phase error 
gives 6(j) = 5x 10^^ radians. Therefore, we expect reflec- 
tion of the outgoing radiation at the outer boundary to 
be insignificant for 30c-l. This is confirmed by the excel- 
lent agreement between runs 30c-l and 30c-2 (the latter 
having even larger outer boundary). 

The second comparison shown in the lower panel of 
Fig. [7] is the phase difference between 30c-l/N6 and 24b- 
1/N4, a shorter 8-orbit evolution started from a sepa- 
rate initial data set (set 24b in Table |T| with a separate 
eccentricity-reduction procedure. The phase agreement 
between these two runs (including an overall time shift 
and phase shift) is better than 0.01 radians for a total 
accumulated phase of ~ 100 radians of the 8-orbit run, 
i.e. better than one part in 10^. Run 24b-l has a similar 
outer boundary location as run 30c-3, and indeed both 
of these runs show similar phase differences from 30c- 1. 

Finally, the third comparison shown in the lower panel 
of Fig. [7] is the phase difference between the two highest 
resolutions of the run 30c-l when a time shift is applied. 
For t > 1000m the agreement is much better than with- 
out the time shift (see upper panel), indicating that the 
dominant error is a small difference in the overall evolu- 
tion time. For the post-Newtonian comparisons we per- 
form in the second part of this paper, waveforms are al- 
ways aligned at specific frequencies by applying time and 
phase shifts. Therefore, the time-shifted phase difference 
as displayed in the lower panel is the most appropriate 
measure of numerical truncation error for these PN com- 
parisons. This difference is less than 0.003 radians after 
t = IOOOto but is larger, about 0.02 radians, at early 
times where the waveforms are noisy because of junk ra- 
diation. 

We now compare the gravitational wave amplitudes 
of different runs in the same manner as we compared the 
gravitational wave phases. Figure[5]presents convergence 
data for the amplitude of the gravitational waves for the 
same runs as shown in Fig. [71 Spatial truncation error 
for the amplitude is less then 0.1 percent for t/m > 1000, 
and earlier than this it is limited by residual noise from 
the junk radiation. Differences (including time shifts) 
between runs of different lengths are shown in the lower 
panel of Fig. [8l These differences are even smaller, but 
because of their small size, they are dominated by noise 
for about the first half of the run. The oscillations appar- 
ent in the comparison to 24b-l are caused by the larger 
orbital eccentricity of 24b- 1 (cf. Tab.[T|. 



F. Extrapolation to infinity 

The quantity of interest to gravitational wave detec- 
tors is the gravitational waveform as seen by an observer 
effectively infinitely far from the source. Our numeri- 
cal simulations, in contrast, cover only a region of finite 
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FIG. 8: Convergence of the gravitational wave amplitude ex- 
tracted at radius R = 77m. This plot corresponds to Fig. (Q, 
except that relative amplitude differences are shown. The 
thin vertical line indicates the time at which mui = 0.1 for 
30C-1/N6. 
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FIG. 9: Difference between areal radius rarcai and coordinate 
radius r of selected extraction surfaces, rarcai remains con- 
stant to within 0.01m during the evolution. The diamond 
indicates AfADM/m of the initial data. 



volume around the source, and our numerical waveforms 
are extracted at a finite radius. Waveforms extracted at 
a finite radius can differ from those extracted at infinity 
because of effects discussed in Section IIIDI these effects 
can lead to phase errors of several tenths of a radian and 
relative amplitude errors of several percent. To avoid 
such errors we extrapolate to infinite extraction radius 
as follows. 

We extract data for '^4 on coordinate spheres of coordi- 
nate radii r/m = 75, 80, 85, ... , 240, as described in Sec- 
tion IIIDI These extracted waveforms are shifted in time 
relative to one another because of the finite light-travel 
time between these extraction surfaces. We correct for 
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FIG. 10: Error of phase extrapolation to infinity for extrapo- 
lation of order n, cf. Eq. (|19|l . Plotted are absolute differences 
between extrapolation with order n and n + 1. Increasing the 
order of the polynomial increases accuracy, but also amplifies 
noise. 



FIG. 11: Error of amplitude extrapolation to infinity for ex- 
trapolation with order n, cf. Eq. (|20|l . Plotted are relative 
amplitude differences between extrapolation with orders n 
and 71+1. The inset is an enlargement for t — r* > 1000m. 



this by shifting each waveform by the tortoise-coordinate 
radius at that extraction point [95| 
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(18) 



Here AfAP M is th e ADM mass of the initial data, and 
''areai = \/ A/Att, where A is the area of the extrac- 
tion sphere. This is not the only possible choice for 
the retarded time — for example, the waveforms could be 
shifted so that the maximaof the amplitude align [41]. It 
has also been suggested (l03j | that the time shift should 
change with the amount of radiated energy — essentially, 
that the factor of A/adm should be replaced by the 
amount of mass interior to the extraction radius at each 
time. We leave investigation of other choices of retarded 
time for future work. 

Figure [9] presents the areal radius during the evolu- 
tion at several typical extraction radii. The areal radius 
of these extraction surfaces is constant to within about 
0.01m, and to the same precision, raroai = r + Madm. 
This relationship is not surprising, because the initial 
data is conformally flat, so that for coordinate spheres 
''areal = T + A/adm + 0{MADM/r). FoT Convenience, we 
simply set raioai = r + Madm in Eq- (|18p. rather than 
explicitly integrating to find the area of each extraction 
sphere. 

After the time shift, each waveform is a function of re- 
tarded time, t — r*. At a given value of retarded time, we 
have a series of data points — one for each extraction ra- 
dius. We flt phase and amplitude of these data separately 



to a polynomial in 1/r, 

c^{t -r\r)^ 0(0) (t - r*) + X: M!^, (19) 



k=l 



rA{t-r\r)=A,,,{t-r*)+j: ^^''^\ (20) 



fc=i 



The leading-order term of each polynomial, as a function 
of retarded time, is then the desired asymptotic wave- 
form: 



</.(*- r* ) = 0(0) (t-r*), 
rA(t-r*) = A(o)(i-r*). 



(21) 
(22) 



We find good convergence of this method as we increase 
the order n of the extrapolating polynomial. Figure 1101 
shows the difference in phase between waveforms extrap- 
olated using successively higher-order polynomials. We 
see a broad improvement in the accuracy of the phase 
with increasing order, but unfortunately, higher order 
extrapolations tend to amplify the noise. Our preferred 
choice is n = 3 extrapolation, resulting in extrapolation 
errors of < 0.003 radians for t — r* > lOOOm. 

Figure [TT] is analogous to Fig. [TOl except that it shows 
relative differences in the extrapolated amplitudes. The 
basic picture agrees with the phase extrapolation: Higher 
order extrapolation reduces the errors, but amplifies 
noise. Our preferred choice n — i gives a relative am- 
plitude error of ^ 0.002 for i — r* > lOOOm, dropping to 
less than 0.001 for t - r* > 2000m. 

Phase and amplitude extrapolation become increas- 
ingly more accurate at late times. The main obstacle 
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FIG. 12: Effect of wave extraction radii on extrapolated 
piiase. Each curve represents the difference from our pre- 
ferred wave extrapolation using r £ [75m, 240m]. The three 
solid curves represent extrapolation from different intervals of 
extraction radii. The curves labeled "240m" and "90m" rep- 
resent differences from waves extracted at these two radii, 
without any extrapolation, for two cases: time and phase 
shifted so that (j) and (p match at mu> = 0.1 (dashed), and 
without these shifts (dotted). 




(t-r*)/m 

FIG. 13: Effect of choice of wave extraction radii on extrap- 
olated amplitude. Each curve represents the (relative) am- 
plitude difference to our preferred wave extrapolation using 
r G [75m, 240m]. The three solid curves represent extrapo- 
lation from different intervals of extraction radii. The curves 
labeled "240m" and "90m" represent differences from waves 
extracted at these two radii, without any extrapolation, for 
two cases: time and phase shifted so that (p and (p match at 
mu} — 0.1 (dashed), and without these shifts (dotted). 



to accuracy seems to be near-zone effects scaling with 
powers of (A/r), where A is the wavelength of the grav- 
itational wave. The wavelength is quite large at the be- 
ginning of the simulation (A w 180m, cf. Fig. O, but 
becomes shorter during the evolution, so that even low- 
order extrapolation is quite accurate at late times. Al- 
ternatively, near-zone effects can be mitigated by using 
data extracted at large values of r. It is precisely be- 
cause of these near-zone effects that we have chosen to 
ignore data extracted at r < 75rn when we extrapolate 
to infinity. 

In Figs. [T2] and [T31 we show the effects of extrapola- 
tion using different ranges of extracted data. Using data 
extracted every 5m in the range r = 50m-90m results in 
noticeable differences early in the run — though it is ad- 
equate later in the run. For ranges at higher radii (e.g. 
[75m, 150m] or [150m, 240m]), the accuracy is not highly 
variable, though we find that noise is increased when us- 
ing data from such a smaller range of extraction radii. 

To estimate the errors generated by not extrapolat- 
ing waveforms to infinity at all. Fig. [T2l contains also the 
phase difference between wave extraction at two finite 
radii (90m and 240m) and our preferred extrapolated 
phase at infinity. The dotted lines show such phase dif- 
ferences when only a time shift by the tortoise-coordinate 
radius of the extraction sphere is applied. The errors are 
dramatic, tenths of radians or more, even very late in the 
run. When matching to post-Newtonian waveforms, we 



are free to add an overall time and phase shift (cf. Sec- 
tion |IVB|. Therefore, the dashed lines in Fig. [T^ show 
phase differences with the same unextrapolated wave- 
forms as shown by the dotted lines, except that a phase 
and time shift has been applied so that the 4> and (p 
agree with those of the extrapolated waveform late in 
the run (where muj — 0.1), where the wavelengths are 
shortest and wave extraction is expected to work best. 
Even with such an adjustment, the gravitational wave 
phase extracted at r = 90m differs by about 0.1 radian 
at t ~ 1000m before coalescence, with this difference 
growing to 0.3 radians at the start of our simulation. 

Figure [13] makes the same comparison for the gravi- 
tational wave amplitude. Wave extraction at r = 90m 
results in relative amplitude errors of up to 8 per cent, 
and of about 2 per cent even in the last 1000m of our sim- 
ulation. We also point out that the errors due to finite ex- 
traction radius decay approximately as the inverse of the 
extraction radius: For waves extracted at r = 240m the 
errors are smaller than for waves extracted at r = 90m 
by about a factor of three, as can be seen in Figs. [T2l 
and [131 for wave extraction at r = 45m, the errors would 
be approximately twice as large as the r = 90m case. 
The errors introduced by using a finite extraction radius 
are significantly larger than our truncation error (even 
at extraction radius 240m). Therefore extrapolation to 
infinity is essential to realize the full accuracy of our sim- 
ulations. 
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G. Estimated time of merger 

Since we have not yet been successful with simulat- 
ing the merger, we do not precisely know when merger 
occurs. However, by comparing the orbital and gravita- 
tional wave frequencies to already published results, we 
can nevertheless estimate the time of merger. 

The simulation presented in Fig. [3] stops at time t = 
3929to when the horizons of the black holes become too 
distorted just before merger. At that point, the proper 
separation between the horizons is ~ 4.0m, and the or- 
bital frequency has reached mfiorbit — 0.125; comparison 
with [3^ suggests this is about 15m before formation of 
a common apparent horizon, i.e. the common horizon 
should form in our simulations at icAH ~ 3945m. 

The waveform extrapolated to infinity ends at t — r* = 
3897to at a gravitational wave frequency of muj w 0.16. 
This places the end of the waveform at about 50m (or 
1.5 cycles) before formation of a common apparent 
horizon^ (judged by comparison with (33|). Thus, we es- 
timate the formation of a common horizon to correspond 
to a retarded time of approximately (i— r*)cAH ~ 3950m. 



III. GENERATION OF POST-NEWTONIAN 
WAVEFORMS 

It is not our intention to review all of post-Newtonian 
(PN) theory, but to summarize the important points that 
go into the construction of the post-Newtonian wave- 
forms that we will compare to our numerical simulation. 
For a complete review of post-Newtonian methods ap- 
plied to inspiralling compact binaries, see the review ar- 
ticle by Blanchetjl04i]. 

The post-Newtonian approximation is a slow-motion, 
weak-field approximation to general relativity with an 
expansion parameter e ^ {v/cf' ~ {Gm/rc^). For a bi- 
nary system of two point masses mi and m2, v is the 
magnitude of the relative velocity, m is the total mass, 
and r is the separation. In order to produce a post- 
Newtonian waveform, it is necessary to solve both the 
post-Newtonian equations of motion describing the bi- 
nary, and the post-Newtonian equations describing the 
generation of gravitational waves. 

Solving the equations of motion yields explicit expres- 
sions for the accelerations of each bod y in terms of the po- 
sitions and vel ocities of the two bo dies [l05l[T06lll07l[l08l 
[T09l . flTol . [ml . im im nil Uil. The two-body equa- 
tions of motion can then be reduced to relative equations 
of motion in the center-of-mass frame in terms of the 



The waveform ends somewhat further from merger than the or- 
bital trajectory, because the artificial boundary is placed initially 
at a radius ~ 15m, and then moves outward somewhat faster 
than the speed of light, thus overtaking the very last part of the 
waveform as it travels to the wave-extraction radii. 



relative position and velocity fTT6'|. The relative accelera- 
tion is currently known through 3.5PN order, where OPN 
order for the equations of motion corresponds to Newto- 
nian gravity. The effects of radiation reaction (due to the 
emission of gravitational waves) enters the relative accel- 
eration starting at 2.5PN order. The relativistic correc- 
tions to the relative acceleration at IPN, 2PN and 3PN 
order (ignoring the radiation reaction terms at 2.5PN and 
3.5PN order) admit a conserved center of mass binding 
energy through 3PN order |ll7| . There is no 2.5PN or 
3.5PN order contribution to the energy. 

Solving the post-Newtonian wave generation problem 
yields expressions for the gravitational waveform hij and 
gravitati onal wave flux C in terms of radiative multipole 
moments [ll8j |. These radiative multipole moments are in 
turn related to the source multipole moments, which can 
be given in terms of the relative position and relative 
velocity of the binary ^11^. For the gravitational wave 
generation problem, PN orders are named with respect 
to the leading order waveform and flux, which are given 
by the quadrupole formalism. Thus, for example, 1.5PN 
order in the wave generation problem represents terms of 
order (v/c)^ beyond quadrupole. Higher order effects en- 
ter both through post-Newtonian corrections to the mass 
quadrupole, as well as effects due to higher multipole mo- 
ments. Starting at 1.5PN order the radiative multipole 
moments include non-linear and non-instantaneous (i.e. 
depend upon the past history of the binary) interactions 
among the source m ultipole m oments (e.g. gravitational 
wave tails) .HE, 120, lUll [HI . 

It was recognized early that simply plugging in the 
orbital evolution predicted by the equations of motion 
into the expressions for the waveform would not gener- 
ate templates accurate enou gh fo r matched filtering in 
detecting gravitational waves [12^. This is because radi- 
ation reaction enters the equations of motion only at the 
2.5PN order; hence computing a waveform to k PN order 
beyond the quadrupole formalism would require 2.5 -I- k 
PN orders in the equations of motion. In order to ob- 
tain as accurate a post-Newtonian waveform as possible 
it is thus necessary to introduce the assumption of an 
adiabatic inspiral of a quasi-circular orbit, as well as the 
assumption of energy balance between the orbital bind- 
ing energy and the energy emitted by the gravitational 
waves. 



A. Adiabatic inspiral of quasi-circular orbits 

The emission of gravitational radiation causes the or- 
bits of an isolated binary system to circularize 50] . Thus 
it is a reasonable assumption to model the orbital evolu- 
tion of the binary as a slow adiabatic inspiral of a quasi- 
circular orbit. With this assumption, post-Newtonian 
expressions for the orbital energy E and gravitational 
energy flux C are curre ntly known through 3.5PN order 
p i il25l . [T26l . [127I IT2i | . These expressions can be given 
in terms of a parameter related to either the harmonic 
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coordinate separation r, or to the orbital frequency Q. 
We choose to use the expressions given in terms of a 
frequency-related parameter 



(23) 



rather than a coordinate-related parameter, because the 
coordinate relationship between the numerical simulation 
and the harmonic coordinates used in post-Newtonian 
approximations is unknown. The orbital energy for an 
equal mass system is given bv |l04| 



E : 



37 1069 2 

X a; 

48 384 

/ 1427365 205 



V 331776 384 



(24) 



and the gra vitat ional wave flux for an equal mass system 
is given by [104 



where 7 = 0.577216 ... is Euler's constant. 



B. Polarization Waveforms 

The gravitational polarization waveforms for a quasi- 
circular orbit in the x — y plane, as measured by an ob- 
server at spherical coordinates (i?, 0, 4>), are given by 

h+ = ^a;|-(l + cos0)cos2($-0) + ---| (26) 



hy = 



c^R' 
2G/i 



|-2cos0sin2($-(^) + ---|, (27) 



where $ is the orbital phase (measured from the x-axis) 
and ^ — 77117712/™ is the reduced mass. The polariz ation 
wav eforms are currently known through 2.5PN order jl29l 
IT30I. 



£ = ^x^ \ 1 - + Anx^/'^ - -^x^ - ^TTX 



2c> 
5G" 



18608019757 355 



59 

567" 
1712 



767 
I2" 



209563200 



856 
105 



In (16a;) 



64 105 
16655 7/9 



6048 



5^2 1- Optimally oriented observer 

For an equal-mass binary the polarization waveforms 
along the 2:-axis (i.e. the optimally oriented observer 
along the normal to the orbital plane) are given by 
(25) pMlT30i 
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where 
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2 , fGm 
3^"te 



(30) 



is a constant frequency scale that depends upon the con- 
stant time scale b entering the gravitation al wave ta il con- 
tribution to the polarization waveforms |l3ll . Il32j . The 
freely-specifiable constant b corresponds to a choice of 
the origin of radiative time T with respect to harmonic 
time t, and enters the relation between the retarded time 
Tr — T — R/c in radiative coordinates (the coordinates 
in which the waveform is given) and the retarded time 
t— r/c in harmonic coordinates (the coordinates in which 



the equations of motion are given) [131l . |13 

2GAfADM 



t- 



■In 



(31) 



Here A/adm is the ADM mass (mass monopole) of the 
binary system. 



2. The (2,2) mode 

When comparing a post-Newtonian waveform with 
data from a physical gravitational wave detector, it is 
necessary to compare waves emitted in a certain direction 
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(0, (f>) with respect to the source. However, comparing 
waveforms between PN and numerical simulations can 
be done in all directions simultaneously by decomposing 
the waveforms in terms of spherical harmonics and then 
comparing different spherical harmonic modes. Since the 
power in each spherical harmonic mode decreases rapidly 
with spherical harmonic index, with the (2, 2) mode dom- 



inating (for an equal-mass non-spinning binary) , it is pos- 
sible to do a very accurate comparison that is valid for 
all angles by using only a few modes. In addition, as 
pointed out by Kidder 62], the dominant (2,2) mode can 
be computed to 3PN order. For an equal-mass binary, 
the (2,2) mode is 
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(32) 



r 



Since the (2,2) mode of the numerical waveforms is less 
noisy than the waveform measured along the z-axis, and 
since we have access to the 3PN amplitude correction of 
the (2,2) mode, we will use the (2,2) waveforms rather 
than the z-axis waveforms for our comparisons between 
NR and PN in Sec. IVII We have verified (for all compar- 
isons using post-Newtonian waveforms of < 2.5PN order 
in amplitude) that our results do not change significantly 
when we use z-axis waveforms instead of (2,2) waveforms. 



these terms are beyond the order to which the orbital 
phase evolution is known (3.5PN order), it can be argued 
that these terms can be ignored. Note that the choices 
of Xq in Eq. (j30p and S in Eq. (j33p are not unique; they 
were made to gather all logarithmic term s into one term, 
as well as to simplify the waveform [133| |. 



D. Energy balance 



C. Absorbing amplitude terms into a redefinition 
of the phase 

The logarithms of the orbital frequency parameter x 
(as well as the constant frequency scale xq) that appear 
in the amplitude expressions (pS)) . P5)) . and ([5^ can be 
absorbed into a redefinition of the phase by introducing 
an auxiliary phase variable 5* = <I> -I- 5. Noting that the 
Inx terms first enter at 1.5 PN order , it is straightforward 
to show that choosing [g^, Il29l Il33j | 



J = -3^^x3/2 In f- 

Xq 



(33) 



where Madm/"i = 1 ~ x/8 + 0{x^) for an equal mass 
system, will eliminate the Inx terms from both the (2,2) 
mode as well as for the polarization waveforms. This 
follows from 



1(2,2) 



= Ai 



-2i<S> 



(l-2i(5-2(52 + 0(a;9/2)) 



and similarly for the polarization waveforms. Further- 
more, since the orbital phase as a function of frequency 
goes as at leading order (see Eq. (|40)) below), the 

In a; terms, which were 1.5PN, 2.5PN, and 3PN order in 
the original amplitude expressions, now appear as phase 
corrections at relative order 4PN, 5PN, and 5.5PN. As 



The second assumption that goes into making as ac- 
curate a post-Newtonian waveform as possible is that of 
energy balance. It is assumed that the energy carried 
away by the emission of gravitational waves is balanced 
by the change in the orbital binding energy of the binary, 



dE 
'dt 



C. 



(34) 



While this is extremely plausible, it has only been con- 
firmed through 1.5 PN order ^134]. 

Given the above expressions for the energy, flux, and 
waveform amplitude, there is still a set of choices that 
must be made in order to produce a post-Newtonian 
waveform that can be compared to our numerical wave- 
form. These include 

1. The PN order through which terms in the orbital 
energy and luminosity arc retained. 

2. The procedure by which the energy balance equa- 
tion is used to obtain x{t) and ^{t). 

3. The PN order through which terms in the waveform 
amplitude are kept. 

4. The treatment of the Inx terms. These terms can 
be included in the amplitude or included in the or- 
bital phase via the auxiliary phase 5* = $ + (5. If 
the latter is chosen, these terms can be retained 
or ignored; ignoring them can be justified because 
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they occur at higher order than all known terms in 
the orbital phase. 

We always expand energy and luminosity to the same 
order, which may be different from the order of the am- 
plitude expansion; both of these expansion orders are 
indicated explicitly in each of our comparisons. We ig- 
nore the ln(a;/a;o) terms in the amplitude by absorbing 
them into the phase and dropping them because of their 
high PN order. In the next section we describe several 
choices for obtaining x{t) and from the energy bal- 
ance equation. 



E. Taylor approximants: Computing 

In this section we describe how to obtain the orbital 
phase as a function of time, using the energy bal- 

ance equation (|34p . Different methods of doing this ex- 
ist; here we follow the naming convention of [5^. These 
methods, and variations of them, are called Taylor ap- 
proximants, and all formally agree to a given PN order 
but differ in how higher-order terms are truncated. We 
discuss four time-domain approximants here, but more 
can be defined. 



1. TaylorTl 

The TaylorTl approximant is obtained by numerically 
integrating the ODEs 



dx 
It 

dt 



{dE/dx) 
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(35) 
(36) 



to produce $(t). The fraction on the right side of Eq ([35]) 
is retained as a ratio of post-Newtonian expansions, and 
is not expanded further before numerical integration. 
This is the approximant used in the NR-PN comparisons 
in [3i,|43. 



2. TaylorT2 

The TaylorT2 approximant is obtained by starting 
with the parametric solution of the energy balance equa- 
tion: 



r° , {dE/dx) 
t{x) = to+ dx^ ' ' 



C 



(37) 



J X 



x^l^c^ (dE/dx) 



The integrand of each expression is re-expanded as a sin- 
gle post-Newtonian expansion in x and truncated at the 



appropriate PN-order; these integrals are then evaluated 
analytically to obtain for an equal-mass binary (s^ . [ssj : 
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3. TaylorTS 

The TaylorT3 approximant is closely related to Tay- 
lorT2. It is obtained by introducing the dimensionless 
time variable 



ito-t), 



(41) 



where v — TOim2 /m2 and r^i/^ = o{e). The TaylorT2 
expression t(x) is inverted to obtain x{t), and truncated 
at the desired PN order. Then x{t) is integrated to ob- 
tain 



$(t) = - / dr 
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This procedure yields for an equal-mass binary [10 
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This is the post-Newtonian approximant used in visual 
comparisons by [s^] and in the NR-PN comparisons in 
[4l| at 3PN order in phase. 
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4. TaylorTi 

In addition to simply numerically integrating the flux- 
energy equation (l35|) . as is done for TaylorTl, one may 
instead re-expand the right side of (|55|) as a single series 
and truncate at the appropriate PN order before doing 
the integration. The phase evolution <i>(t) can thus be 
obtained by numerically integrating the ODEs 



dx 
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This approximant was not considered in |54j . however 
for consistency with their notation, we call it TaylorT4. 
TaylorT4 is the primary approximant used in the NR-PN 
comparisons in [38, 39], and one of the several approx- 
imants considered in the NR-PN comparisons in [ssj . 
Ref.t32i] pointed out that TaylorT4 at 3.5PN order in 
phase is close to TaylorT3 at 3PN order in phase, and 
therefore should give similar agreement with numerical 
results. 



IV. PN-NR COMPARISON PROCEDURE 
A. What to compare? 

There are many ways to compare numerical relativ- 
ity and post-Newtonian results. For example, the post- 
Newtonian orbital phase ^{t) could be compared with 
the coordinate phase of the black hole trajectories. How- 
ever, this and many other comparisons are difficult to 
make in a coordinate-independent manner without ex- 
pending significant effort to understand the relationship 
between the gauge choices used in post- Newtonian theory 
and in the NR simulations. Therefore, in order to obtain 
the most meaningful comparison possible, we attempt to 
minimize gauge effects by comparing gravitational wave- 
forms as seen by an observer at infinity. The waveform 
quantity most easily obtained from the numerical rela- 
tivity code is the Newman-Penrose quantity ^'4, and we 
will compare its (2,2) component [cf. Eq. P^ ]. split 
into phase and amplitude A according to Eq. ([M]) and 
extrapolated to infinite extraction radius. 

The post- Newtonian formulae in Section |TTT] yield the 
metric perturbation components and hx, which — for 
a gravitational wave at infinity — are related to 5*4 by 



*4(i) 



We numerically differentiate the post-Newtonian expres- 
sions for h+(t) and hx{t) twice before computing ampli- 
tude and phase using Eq. (fT4| . Note that (j)(t) will differ 
slightly from the phase computed from the metric per- 
turbation directly, as tan~^(ft,x/ft.+), because both the 
amplitude and phase of the metric perturbation are time 
dependent. For the same reason, (j){t) is not precisely 
equal to twice the orbital phase. 

As in Ref. [4l[, we compare ^'4 rather than /i+.x 
to avoid difficulties arising with fixing the integration 
constants when integrating the numerically obtained ^'4 
(see [13] for more details). Both 5*4 and /i+.x contain 
the same information, so differences between both proce- 
dures should be minimal. 



B. Matching procedure 

Each of the post-Newtonian waveforms has an arbi- 
trary time offset to and an arbitrary phase offset 0o- 
These constants can be thought of as representing the 
absolute time of merger and the orientation of the bi- 
nary at merger, and we are free to adjust them in order 
to match NR and PN waveforms. Following jH, |4l[, we 
choose these constants by demanding that the PN and 
NR gravitational wave phase and gravitational wave fre- 
quency agree at some fiducial frequency LUm ■ Specifically, 
we proceed as follows: We start with a NR waveform 
^'^^ (i) and an unshifted PN waveform 5'f^'(i) that has 
an arbitrary time and phase shift. After selecting the 
matching frequency Wm, we can find (to essentially un- 
limited accuracy) the time tc such that the derivative of 
the PN phase satisfies 0pN'(ic) = i-^rm where 4>p]S!'{t) is 
the phase associated with ^'4^ {t). Similarly, we find the 
time tm such that 0NR(im) = '^m- The time tm cannot be 
found to unlimited accuracy, and the uncertainty in tm is 
due mainly to residual eccentricity of the NR waveform, 
as discussed in Section IV El Once we have t„i and tc, we 
leave the NR waveform untouched, but we construct a 
new, shifted, PN waveform 



i(0PN' (tc)-0NR(tm)) 



(48) 



(h+it) - ihxit)) ^ 
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The phase of this new PN waveform is therefore 

(/'Pn(<) = (pPN'it + tc~ tm) - 0PN'(tc) + 0NR(im), (49) 

which satisfies 0pN(im) = 0NR(im) and ^pN(tm) = i^m 
as desired. All our comparisons are then made using the 
new shifted waveform \E'4^(i) rather than the unshifted 
waveform (t). 



C. Choice of Masses 

The post-Newtonian expressions as written in Sec- 
tion mil involve the total mass m, which corresponds to 
the the sum of the bare masses of the point particles in 
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post-Newtonian theory. When comparing PN to NR, the 
question then arises as to which of the many definitions of 
the mass of a numerically-generated binary black hole so- 
lution should correspond to the post-Newtonian param- 
eter m. For non-spinning black holes at very large sepa- 
ration, m reduces to the sum of the irreducible masses of 
the two holes. Neglecting tidal heating, the irreducible 
masses should be conserved during the inspiral, so that 
we identify m with the sum of the irreducible masses of 
the initial data 30c. As discussed in Sec. |V] the black 
hole spins are sufficiently small so that there is no dis- 
cernible difference between irreducible mass of the black 
holes and the Christodoulou mass, Eq. ([3]). Of course, 
the latter would be more appropriate for spinning black 
holes. 



V. ESTIMATION OF UNCERTAINTIES 

To make precise statements about agreement or dis- 
agreement between numerical and post-Newtonian wave- 
forms, it is essential to know the size of the uncertainties 
in this comparison. When discussing these uncertainties, 
we must strive to include all effects that may cause our 
numerical waveform to differ from the post-Newtonian 
waveforms we compare to. For instance, in addition to 
considering effects such as numerical truncation error, 
we also account for the fact that NR and PN waveforms 
correspond to slightly different physical scenarios: The 
PN waveforms have identically zero spin and eccentric- 
ity, whereas the numerical simulations have some small 
residual spin and eccentricity. Table Hill lists all effects we 
have considered; we discuss these in detail below start- 
ing in Sec. IV Al All uncertainties are quoted in terms of 
phase and amplitude differences, and apply to waveform 
comparisons via matching at a fixed ujm according to the 
procedure in Sec. IIVBI 

Most of the effects responsible for our uncertainties are 
time dependent, so that it is difficult to arrive at a single 
number describing each effect. For simplicity, the error 
bounds in Table Hill ignore the junk-radiation noise that 
occurs in the numerical waveform for t — r* < 1000m. 
The extent to which this noise affects the PN-NR com- 
parisons presented below in Sections IVI Al and IVI Bl will 
be evident from the noise in the graphs in these sections. 
Note that all four matching frequencies ujm occur after 
the noise disappears at i — r* ~ 1000m. Furthermore, 
the post-Newtonian waveforms end at different times de- 
pending on the PN order and on which particular post- 
Newtonian approximant is used. Therefore, in order to 
produce a single number for each effect listed in Table HIH 
we consider only the part of the waveform prior to some 
cutoff time, which we choose to be the time at which the 
numerical waveform reaches gravitational wave frequency 
muj = 0.1. 



TABLE III: Summary of uncertainties in the comparison be- 
tween numerical relativity and post-Newtonian expansions. 
Quoted error estimates ignore the junk-radiation noise at 
t < 1000m and apply to times before the numerical waveform 
reaches gravitational wave frequency mu = 0.1. Uncertainties 
apply to waveform comparisons via matching at a fixed Um 
according to the procedure in Sec. IIV Bl and represent the 
maximum values for all four different matching frequencies 
ujm that we consider, unless noted otherwise. 



Effect c 


(radians) 5 A/ A 


Numerical truncation error 


0.003 


0.001 


Finite outer boundary 


0.005 


0.002 


Extrapolation r oo 


0.005 


0.002 


Wave extraction at rarcai=const? 


0.002 


10"* 


Drift of mass m 


0.002 


10-* 


Coordinate time = proper time? 


0.002 


10"* 


Lapse spherically symmetric? 


0.01 


4 X lO"* 


residual eccentricity 


0.02" 


0.004 


residual spins 


0.03 


0.001 


root-mean-square sum 


0.04" 


0.005 



"For the case of matching at mUm = 0.04, the phase uncertainty 
due to residual eccentricity increases to 0.05 radians, thus increas- 
ing the root-mean-square sum to 0.06 radians. 



A. Errors in numerical approximations 

The first three error sources listed in Table HIJ have al- 
ready been discussed in detail in Section [TTj We estimate 
numerical truncation error using the difference between 
the two highest resolution runs after the waveforms have 
been shifted to agree at some matching frequency ujm- 
For mujm —0.1 this difference is shown as the curves la- 
beled '30c-l/N5' in the lower panels of Figs. [7] and [51 and 
corresponds to a phase difference of 0.003 radians and a 
relative amplitude difference of 0.001. For other values 
of u!m the differences are similar. The effect of the outer 
boundary is estimated by the difference between the runs 
30c-l/N6 and 30c-2/N6, which for muim — 0.1 is shown 
as the curves labeled '30c-2/N6' in the lower panels of 
Figs. [71 and [51 and amount to phase differences of 0.005 
radians and relative amplitude differences of 0.002. Er- 
rors associated with extrapolation to infinity have been 
discussed in detail in Figs.[T0land[T2l Specifically, Fig.lTOl 
shows that increasing the extrapolation order between 3 
and 4 changes the extrapolated phase by less than 0.005 
radians, and Fig. [T^confirms that the extrapolated result 
is robust under changes of extraction radii. 

B. Constancy of extraction radii 

If the physical locations of the coordinate-stationary 
extraction radii happen to change during the evolution, 
then the extracted gravitational waves will accrue a tim- 
ing error equal to the light-travel time between the orig- 
inal location and the final location. From Fig. [SI we see 



20 



that the drift in areal radius is less than 0.02m, resulting 
in a time uncertainty of St — 0.02to. This time uncer- 
tainty translates into a phase uncertainty via 



muj X (5t/m) 



(50) 



which yields 5(j) ~ 0.002, when mw = 0.1 (the value at 
the end of the PN comparison) was used. 

To estimate the effect of this time uncertainty on the 
amplitude, we first note that to lowest order in the post- 
Newtonian parameter x (defined in Eq. (|23p ). the wave 
amplitude of ^4 scales like x'^. Also, from Eq. ([l5|) . we 
have dx/dt — 16/(5m)x^. Therefore, 



SA 
A 



dlnAdx 
dx dt 



5t^^(^^/2)8/3^, (51) 

5 TO 



where we have used the fact that the gravitational wave 
frequency uj is approximately twice the orbital frequency. 
For a time uncertainty 5t = 0.02to, Eq. (|?T|) gives 
5A/A « 10-4 for muj = 0.1. 



C. Constancy of mass 

Our comparisons with post-Newtonian formulae as- 
sume a constant post-Newtonian mass parameter to, 
which we set equal to the total irreducible mass of the 
black holes in the numerical simulation. If the total 
mass of the numerical simulation is not constant, this 
will lead to errors in the comparison. For example, 
changes vnt/m caused by a changing mass will lead to 
phase differences. Figure |4] demonstrates that the ir- 
reducible mass is conserved to a fractional accuracy of 
about 5m /m « 5 x 10~®. 

This change in irreducible mass could be caused by nu- 
merical errors, or by a physical increase of the mass of 
each black hole through tidal heating. For our simula- 
tions, m{t) decreases during the run (this is not apparent 
from Fig. [5] which plots absolute values) , thus contra- 
dicting the second law of black hole thermodynamics. 
Moreover, the increase in m{t) through tidal heating is 
much sma ller than the observed variations in TO(t) (see, 
e.g. [l35l |). Therefore, the variations in m{t) are numer- 
ical errors, and we need to bound the influence of these 
errors on the comparison to post-Newtonian expansions. 

Over an evolution time of t/m = 4000, the observed 
mass uncertainty of 6m/ m « 5 x 10~^ results in an 
uncertainty in the overall time interval of 5{t/m) = 
(t/m) X {5m/ m) w 0.02. This time uncertainty translates 
into a phase uncertainty of Scf) k, 0.002, using Eq. ([50)1 
for muj — 0.1. Note that the effect of the black-hole 
spins on the mass is negligible relative to the numerical 
drift of 5 X 10~^. This is because the spins of the holes 
are bounded by S/Min < 2 x 10^^ and the spin enters 
quadratically into the Christodoulou formula The 
error in the gravitational wave amplitude caused by time 



uncertainties due to varying mass is 5 A/ A 
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FIG. 14: Asymptotic behavior of the lapse at large radii for 
times t/m = 0, 1900, 3800. The top figure displays the an- 
gular average of the lapse as a function of radius ai t — 
0, 1900m, 7800m. The bottom figure shows the dominant 
higher multipole moments of the lapse. Both horizontal axes 
are spaced in 1/r. 



affect the amplitude not only via a time offset, but also 
because the amplitude is proportional to (ci;to/2)*/'^ (to 
lowest PN order). However, this additional error is very 
smah, 5A/A « {8/3)5m/m « 10"^. 



D. Time coordinate ambiguity 

We now turn to two possible sources of error that have 
not yet been discussed, both of which are related to am- 
biguity in the time coordinate. The basic issue is that 
the time variable t in post-Newtonian expansions cor- 
responds to proper time in the asymptotically flat re- 
gion, but the time t in numerical simulations is coordi- 
nate time. These two quantities agree only if the lapse 
function N approaches unity at large distances. To verify 
this, we decompose N in spherical harmonics centered on 
the center of mass of the system. 



(52) 



/=0 m=-/ 



ing Eq. (HTD for 



0.1. An error in the mass will 



The angular average of the lapse function, N{r) = 
^/^nNooir) should then approach unity for r ^ 00, and 
all other modes Nim{r) should decay to zero. The top 
panel of Fig. [T4l plots N{r) — 1 vs. m/r for three differ- 
ent evolution times. Fitting N{r) — 1 for r > IOOto to 
a polynomial in m/r gives a y-intercept of < 5 x 10~^ 
for all three times, and for polynomial orders of two 
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through five. Therefore, the coordinate time of the evo- 
lution agrees with proper time at infinity to better than 
5t/m = t/m X 5 X 10~^ w 0.02, which induces a phase 
error of at most Scj) ~ 0.002 and an amphtude error of 
SA/A « 10-4 [cf. Eqs. ^ and dST])]. 

The second source of error related to the lapse is shown 
in the lower panel of Fig. [Ml which presents the three 
dominant higher order moments Nim (r) . All these modes 
decay to zero as r — > oo, except, perhaps, the real part 
of the N22 mode at t/m = 3800. This mode seems to 
approach a value of about 5 x 10"^. At t = 1900m, 
this mode still decays nicely to zero, hence the maximum 
time uncertainty introduced by this effect at late times 
is St = 1900m X 5 X 10^^ « 0.1m, resulting in a potential 
phase uncertainty of 6(j) ~ 0.01 and a potential amplitude 
uncertainty of 6A/A ^ 4 x 10~^. 



E. Eccentricity 

We estimated the eccentricity during the numerical 
simulation with several of the methods described in 
[3^ . liol . [6l| . and have found consistently e 6 x 10"^. 
This eccentricity can affect our comparison to a post- 
Newtonian waveform of a quasi-circular (i.e. zero eccen- 
tricity) inspiral in three ways. 



1. Change in rate of inspiral 

The first effect arises because an eccentric binary has a 
different inspiral rate than a non-eccentric binary; physi- 
cally, this can be understood by noting that the gravita- 
tional flux and orbital energy depend upon the eccentric- 
ity, and therefore modify the rate at which the orbita l fre- 
quency evolves assuming energy balance. Reference |l36j | 
has derived the first-order correction in the phase of the 
gravitational wave due to this effect. Converting their 
result to our notation and restricting to the equal mass 
case yields 
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5Gm 



(dx/dt) IGc^x'^ 



157 2 f 



19/6 



(53) 



where is the initial eccentricity and Xi is the initial 
value of the orbital frequency parameter. Substituting 
this into Eq. ^ yields 
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17/3 



(54) 



Using Ci = 6 X 10"^ and integrating over the frequency 
range from the start of our simulation to the matching 
frequency of mw = 0.1 yields a phase shift of^^ — 2xl0~^, 
which is dwarfed by many other error sources, such as the 
uncertainty in the numerical mass m, cf. Sec. IV CI 



2. Uncertainty in matching time 

The second way in which eccentricity affects our com- 
parison is by introducing errors in our procedure for 
matching the PN and NR waveforms. Recall that our 
matching procedure involves determining the time tm at 
which the gravitational wave frequency lo takes a cer- 
tain value mojm] eccentricity modulates the instanta- 
neous gravitational wave frequency uj{t) via 



uj{t) = uo{t) [1 + 2e cos(ri^i)] , 



(55) 



where uj(t) represents the averaged "non-eccentric" evo- 
lution of the gravitational wave frequency, and fi^ is the 
frequency of radial oscillations, which is approximately 
equal to the orbital frequency. We see that lo can differ 
from (D by as much as 2eu) ~ 2euj. This could induce an 
error in the determination of tm by as much as 



\St I = 1^ ^ ^ 

LU Co 



(56) 



We can simplify this expression by using Eq. (j45p to low- 
est order, and by noting that the gravitational wave fre- 
quency is approximately twice the orbital frequency. We 
find 



\5tr 



< e- 



bm fmuj 



12 V 2 



-8/3 



(57) 



This uncertainty is largest at small frequencies, because 
the frequency changes much more slowly. For mw = 0.04, 
we find \5tm\ ^ 0.9m, and for muj — 0.1, we find \5tm\ ^ 
0.1m. 

To determine how uncertainties in i,„ translate into 
phase differences, recall that in the matching procedure 
described in Sect ion HVBI i„ enters into the phase of the 
shifted PN waveform according to Eq. Therefore 
the phase difference that we compute between the PN 
and NR waveforms is 

A(/)(t) = 0PN(i) -(/-NRit) 

= 4'm'{t + tc-tm) - 0Nr(O + 0NR(*m) " 4>Vw{tc)- 

(58) 

Then the error in A(/) is found by Taylor expanding 
Eq. jSll): 



54> = (5(A^(i)) = {4>Y>N'{t + tc- tm) - 0NR(im)) Stm. 

= (0PN(t) - ^m^ 6tm- (59) 

Our simulations (and therefore the comparisons to 
post-Newtonian theory) start at moj « 0.033, so that the 
maximal error 6(f> within our comparison at times before 
the matching frequency will be 



ll^^bcforol < |0.033-Wm| \Stjn\ 



(60) 



Combining Eqs. (pO)) and and using e w 6 x 10 ^, 
we find that (5(/)bcforc < 0.01 radians for all four of our 
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matching frequencies mu>m = 0.04,0.05,0.063,0.1. The 
maximum error Sip within our comparison at times after 
the matching frequency is 



l^^aftcrl < |0.1 



\Str. 



(61) 



because we end our comparisons to post-Newtonian the- 
ory at mui = 0.1. Eq. (pTj) evaluates to 0.05 radians for 
mujm = 0.04, and is less than about 0.02 radians for the 
three higher matching frequencies. 

The error in the gravitational wave amplitude caused 
by an error in tm can be estimated by Eq. (j5ip . A conser- 
vative estimate using St — 0.9m still gives a small error, 
SA/A « 0.004. 

Note that the bounds on 50bcforo and (50aftcr are pro- 
portional to the eccentricity of the numerical simulation. 
Even with eccentricity as low as 6 x 10~^, this effect is one 
of our largest sources or error for the PN-NR comparison, 
(cf. Table |TTT|). This is the reason why the simpler eccen- 
tricity removal procedure of Husa et al. [6l| (resulting in 
e = 0.0016) is not adequate for our purposes. 



pute spin-orbit coupling up to 2.5 post-Newtonian order, 
and find that the orbital phase, Eq. (|40p . acquires the 
following spin contributions 
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(63) 



where Xi = ' L/m^ is the projection of the dimension- 
less spin of the i-th hole onto the orbital angular momen- 
tum. For equal-mass binaries with spins xi = X2 = Xi 
this reduces to 
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(64) 



3. Periodic modulation of phase and amplitude 

The third effect of orbital eccentricity is a periodic 
modulation of the gravitational wave phase and ampli- 
tude. If we assume that u}{t) varies on much longer time 
scales than l/fi^ (which is true at large separation) then 
time integration of Eq. ([55]) yields 



2e— sin(rir^)- 



(62) 



Because fir ~ ^ ~ <^/2, we therefore find that the 
gravitational wave phase consists of the sum of the de- 
sired "circular" phase, (j){f), plus an oscillatory compo- 
nent with amplitude 4e « 2 x 10"''. This oscillatory 
component, however, is much smaller than other uncer- 
tainties of the comparison, for instance the uncertainty 
in determination of tm- 

Residual eccentricity will also cause a modulation of 
the gravitational wave amplitude in a manner similar 
to that of the phase. This is because eccentricity ex- 
plicitly ent ers t he post-Newtonian amplitude formula at 
OPN order [l37j | . This term is proportional to e, and since 
e < 6 X 10~^ its contribution to the amplitude error is 
small compared to the effect due to uncertainty in tm- 

While oscillations in phase and amplitude due to ec- 
centricity are tiny and dwarfed by other uncertainties in 
the PN-NR comparison, their characteristic oscillatory 
behavior makes them nevertheless visible on some of the 
graphs we present below, for instance, both panels of 

Fig.iini 



F. Spin 

We now turn our attention to effects of the small resid- 
ual spins of the black holes. References |l38l . Il39j | com- 



Our comparisons to post-Newtonian theory are per- 
formed over the orbital frequency range of 0.0167 < 
mVl < 0.05, corresponding to 0.065 < a; < 0.136. The 
gravitational wave phase is approximately twice the or- 
bital phase, so that the spin-orbit coupling contributes 



= 2[$s(0.065) - $s(0.136) 



-64 X 



(65) 



to the gravitational wave phase. In Sec. lIIBI we estimated 
|S|/M^j. < 5 X lO""*, where Mi„ is the irreducible mass of 
either black hole. Because x < \S\/M?rr ~ 5 x 10"'^, the 
residual black hole spins contribute less than 0.03 radians 
to the overall gravitational wave phase. 

We now turn to errors in the amplitude comparison 
caused by residual spin. From Eq. (j64p we can compute 
the error in orbital frequency as 



sn = 
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(66) 



where we have used Eq. pS]) . Because the amplitude of 
\I'4 scales like fl^/"^, we arrive at 



SA_8Sn_ 128 /235 _^ 270625\ 
~A"3ll~^^ Tb \ 96^ ^ 16128 J 



(67) 



which for mujm = 0.1 (i.e. x = 0.136) gives SA/A = 
2.Qx ^ 1.0 X lO"'''. 

Spi n-orbit te rms also contribute directly to the ampli- 
tude [140l Il4l| . The leading order contribution (for an 
equal-mass binary with equal spins) contributes a term 
SA/A - (4/3)xa;^/^, which is the same order of magni- 
tude as the previous error, 10~^. 
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FIG. 15: Comparison of numerical simulation with TaylorTl 3.5/2.5 waveforms. Left: Difference in gravitational wave phase. 
Right: Relative amplitude difference. Plotted are comparisons for four values of uim- The filled diamond on each curve shows 
the point at which 4> = '-^m- The insets show enlargements for small differences and early times. Also shown is the difference 
between the numerical and restricted (i.e. 3.5PN phase, OPN amplitude) Taylor Tl for mujm ~ 0.1. 



VI. RESULTS 

A. Comparison with individual post-Newtonian 
approximants 

We compare our simulations with four different post- 
Newtonian approximants: the TaylorTl, TaylorT2, Tay- 
lorTS, and TaylorT4 waveforms. These four wave- 
forms agree with each other up to their respective post- 
Newtonian expansion orders, but they differ in the way 
that the uncontrolled higher order terms enter. We start 
with the comparison to TaylorTl. 



1. TaylorTl (3.5PN phase, 2. 5PN amplitude) 

Figure [15] compares the numerical simulation to Tay- 
lorTl 3.5/2.5 waveforms (i.e. expansion order 3.5PN in 
phase and 2.5PN in amplitude, the highest expansion or- 
ders currently available for generic direction, cf. IIIIB[) . 
The left panel shows the phase difference, where we find 
differences of more than a radian for all four matching 
frequencies we consider: LUm = 0.04, 0.05, 0.063, and 
0.01. 

For our largest matching frequency, mujrn = 0.1, the 
phase differences are small toward the end of the run by 
construction. Nevertheless, a phase difference of more 
than 0.5 radians builds up in the ~ 1.5 cycles after the 
matching point before the TaylorTl template generation 
fails. Recall that majm — 0.1 occurs about 2.2 gravita- 
tional wave cycles before our simulations fail, which is 
still about 1.5 cycles before merger. However, the largest 
phase disagreement for mujrn = 0.1 builds up at early 
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FIG. 16: Numerical and TaylorTl 3.5/2.5 waveforms. The 
PN waveform is matched to the numerical one at muim ~ 0.04, 
indicated by the small circle. The lower panel shows a detailed 
view of the last 10 gravitational wave cycles. 



times, reaching 1.5 radians at the beginning of our simu- 
lation, about 28 cycles before the matching 30 cycles 
before the end of the simulation), and still showing no 
sign of flattening even at the start of our simulation. 

To achieve phase coherence with the early inspiral 
waveform, it is therefore necessary to match earlier than 
mojm =0.1. The left panel of Fig. [15] clearly shows that 
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FIG. 17: Comparison of numerical simulation with TaylorT2 3.5/2.5 waveforms. Left: Difference in gravitational wave phase. 
Right: Relative amplitude difference. Plotted are comparisons for four values of Um- The filled diamond on each curve shows 
the point at which (f) = lj,„. The insets show enlargements for small differences and early times. 



phase differences at earlier times become smaller when 
the matching point itself is moved to earlier time. For 
instance, = 0.063 (about eight gravitational wave 

cycles before the end of our simulation), results in phase 
differences less than 0.5 radians during the 22 earlier 
cycles of our evolution. However, the phase difference 
0PN — 0NR does not level off at early times within the 
length of our simulation, so it seems quite possible that 
the phase difference may grow to a full radian or more 
if the numerical simulations could cover many more cy- 
cles. We thus estimate that for TaylorTl, to achieve 
1-radian phase coherence with the early inspiral may re- 
quire matching more than 10 cycles before merger. To 
achieve more stringent error bounds in phase coherence 
will require matching even earlier: for instance it appears 
one needs to use ma;,„ = 0.04 (about 20 cycles before the 
end of our simulation) for a phase error of less than 0.1 
radians. 

While matching at small cj,„ yields good phase coher- 
ence early in the run, it produces much larger phase 
differences late in the run. For example, matching at 
mUm = 0.04 results in a phase difference of almost 2 
radians at frequency muj = 0.1. This rather dramatic 
disagreement is illustrated in Fig. [111 which plots both 
the numerical and the TaylorTl waveform, matched at 
muJm = 0.04. 

The left panel of Fig. [15] also includes a comparison to 
the so-called restricted TaylorTl template, where only 
the leading order amplitude terms are used (i.e. OPN 
in amplitude). The reason that higher-order amplitude 
terms affect the phase differences at all is because we 
are plotting gravitational-wave phase, not orbital phase. 
However, we see that the effect of these higher-order am- 
plitude terms on the phase difference is small. 

We now turn our attention to comparing the ampli- 



tudes of the post-Newtonian and numerical waveforms. 
The right panel of Fig. [TSl shows relative amplitude differ- 
ences between TaylorTl 3.5/2.5 and the numerical wave- 
forms. At early times, the amplitudes agree to within 2 
or 3 per cent, the agreement being somewhat better when 
the matching is performed at early times. At late times, 
the amplitudes disagree dramatically; a large fraction 
of this disagreement lies probably in the fact the post- 
Newtonian point of merger (i.e. the point at which the 
amplitude diverges) occurs at a different time than the 
numerical point of merger. We also plot the amplitude 
of the restricted TaylorTl template. The disagreement 
between restricted TaylorTl and the numerical result is 
much larger, about 5 per cent. 



Hannam et al. [4]J performed a similar comparison, 
matching their waveforms with a restricted TaylorTl 
waveform (i.e. 3.5/0.0) generated using the LIGO Al- 
gorithm Library (LAL) [l42j |. The phase difference they 
observe for waveforms matched at mto = 0.1 is consis- 
tent with our results within numerical errors. When 
matching TaylorTl 3.5/0.0 early in their simulation (at 
muj = 0.0455), however, Hannam et al. find a cumulative 
phase difference of 0.6 radians at muj = 0.1. From Fig.fTSl 
we find a quite different value of 1.5 radians for our sim- 
ulation. This disagreement might be caused by the use 
of the finite extraction radius R — 90m for the gravita- 
tional wave phase in Hannam et. al.: Figure [TU] shows 
that extracting at a finite radius leads to a systematic 
phase error, which will induce a systematic error in de- 
termination of the matching time of Hannam et al. This 
error is amplified by the increasing gravitational wave 
frequency toward merger. 
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FIG. 18: Comparison of numerical simulation with TaylorTS 3.5/2.5 waveforms. Left: Difference in gravitational wave phase. 
Right: Relative amplitude difference. Plotted are comparisons for three values ujrn- The filled diamond on each curve shows 
the point at which (f> — ujm- The lines end when the frequency of the TaylorTS waveform reaches its maximum, which happens 
before mui = 0.1, so that the matching frequency mujm = 0.1 is absent. The left plot also contains TaylorTS 3.0/S.O, matched 
at muj,n — 0.1. The insets show enlargements for small differences. 



2. TaylorT2 (3.5PN phase, 2. 5PN amplitude) 

Figure [T7l presents the comparison between the numer- 
ical waveform and the TaylorT2 approximant. The over- 
all trends are very similar to the TaylorTl comparison 
of Fig. I15i however, the phase differences are smaller by 
about a factor of 2 when matching at mujm = 0.1, and 
smaller by a factor of 3 to 4 when matching earlier. To 
our knowledge TaylorT2 has never been compared to a 
numerical simulation; we include it here mainly for com- 
pleteness. 



3. TaylorTS (3.5PN and 3.0PN phase, 2.5PN amplitude) 

Figure [TS] is the same as Fig US] except it compares 
numerical simulations to the TaylorTS family of wave- 
forms. Two differences between TaylorTl and TaylorTS 
are readily apparent from comparing these two figures. 
The first is that we do not match TaylorTS S.5/2.5 wave- 
forms at muim = 0.1. This is because the frequency of 
TaylorTS S.5/2.5 waveforms reaches a maximum shortly 
before the formal coalescence time of the post-Newtonian 
template, and then decreases. The maximal frequency is 
less than 0.1, so that matching at niLOm = 0.1 is not pos- 
sible. For this reason, we have also shown in Fig. [TH] a 
comparison with a TaylorTS S.O/S.O waveform matched 
at mujm = 0.1. The other major difference between 
the TaylorTS S.5/2.5 and TaylorTl S.5/2.5 comparison is 
that the phase difference, 0pn — '?^'nr, has a different sign. 
While TaylorTl S.5/2.5 spirals in more rapidly than the 
numerical simulation, TaylorTS S.5/2.5 lags behind. In- 



terestingly, the phase differences from the numerical sim- 
ulation for both TaylorTl S.5/2.5 and TaylorTS S.5/2.5 
are of about equal magnitude (but opposite sign). The 
TaylorTS S.O/S.O comparison matched at mujm = 0.1 has 
smaller phase differences than does the TaylorTS 3.5/2.5 
comparison, but the slope of the S.O/S.O curve in Fig. [T51 
is nonzero at early times, so it appears that Taylor TS 
S.O/S.O will accumulate significant phase differences at 
even earlier times, prior to the start of our simulation. 
In Fig.[221it can be seen that matching TaylorTS S.O/S.O 
at mojm — 0.04 leads to a good match early, but leads to 
a phase difference of 0.6 radians by muj = 0.1. 

Hannam et al. [4l| match a TaylorTS S. 0/0.0 wave- 
form at mujm = 0.1 and mujm — 0.0455. Matching 
at mujm — 0.1 again gives phase differences consistent 
with our results within numerical errors. Matching at 
muJm — 0.0455, Hannam et al. find a phase difference 
of 0.9 radians, while we find a smaller value of 0.5 ra- 
dians. Again, this difference could be due to the finite 
extraction radius used by Hannam et al. 

4. TaylorTi (3.5PN phase, 2. 5PN amplitude) 

Figure [in] is the same as Figs . [T5I and [T5I except it com- 
pares numerical simulations to the TaylorT4 FN wave- 
forms. The agreement between TaylorT4 waveforms and 
the numerical results is astonishingly good, far better 
than the agreement between NR and either TaylorTl or 
TaylorTS. The gravitational wave phase difference lies 
within our error bounds for the entire comparison region 
moj < 0.1, agreeing to 0.05 radians or better over 29 
of SO gravitational wave cycles. Ref. jsl] found agree- 
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FIG. 19: Comparison of numerical simulation with TaylorT4 3.5/2.5 waveforms. Left: Difference in gravitational wave 
phase. Right: Relative amplitude difference. Plotted are comparisons for four values of curn- The filled diamond on each curve 
shows the point at which (p — uj,„. The left plot also includes two phase comparisons with expansions of different PN order in 
amplitude, as labeled, for mcjrn = 0.1. 



ment between TaylorT4 and their numerical simulation 
to the level of their numerical accuracy 2 radians), 
agreeing to roughly 0.5 radians in the gravitational fre- 
quency range of 0.054 < muj < 0.1. Ref. 35] found that 
NR agrees better with TaylorT4 than with TaylorTl, but 
the larger systematic and numerical errors of the numer- 
ical waveforms used in these studies did not allow them 
to see the surprising degree to which NR and TaylorT4 
agree. The gravitational wave amplitude of TaylorT4 
agrees with the NR waveform to about 1-2 percent at 
early times, and 8 percent at late times. In Fig. [201 we 
plot the NR and TaylorT4 waveforms; the two waveforms 
are visually indistinguishable on the plot, except for small 
amplitude differences in the final cycles. 

On the left panel of Fig. [19] we also show phase com- 
parisons using PN waveforms computed to 3.5PN order 
in phase but to OPN and 3.0PN orders in amplitude, 
for the case mujm = 0.1. The PN order of the ampli- 
tude expansion affects the phase comparison because we 
are plotting differences in gravitational-wave phase and 
not orbital phase. The differences between using OPN, 
2.5PN, and 3. OPN amplitude expansions are evident on 
the scale of the graph, but because these differences are 
smaller than our estimated uncertainties (see Table IIIII) , 
we cannot reliably conclude which of these most closely 
agrees with the true waveform. 

Figure [21] presents amplitude differences between NR 
and TaylorT4 as the post-Newtonian order of the am- 
plitude expansion is varied, but the phase expansion re- 
mains at 3.5PN. The 2.5PN amplitude curve was already 
included in the right panel of Fig. [19] We see clearly 
that higher order amplitude corrections generally result 
in smaller differences. The 3PN amplitude correction 
to the (2, 2) mode recently derived by Kidder [6^ im- 
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FIG. 20: Numerical and TaylorT4 3.5/3.0 waveforms. The 
PN waveform is matched to the numerical one at muim ~ 0.04, 
indicated by the small circle. The lower panel shows a detailed 
view of the end of the waveform. 



proves agreement dramatically over the widely known 
2.5PN amplitude formulae. Unfortunately, the 3PN am- 
plitude correction to the entire waveform, including all 
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FIG. 21: TaylorT4 amplitude comparison for different PN 
orders. Shown is the relative difference in gravitational wave 
amplitude between TaylorT4 and numerical 122 waveforms as 
a function of time. Matching is performed at mujm ~ 0.04. 
All curves use 3.5PN order in phase but different PN orders 
(as labeled) in the amplitude expansion. 
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FIG. 22: Phase comparison for different PN approximants 
at different PN orders, matched at muirn ~ 0.04. Shown is 
the difference in gravitational wave phase between each post- 
Newtonian approximant and the numerical 122 waveforms as 
a function of time. The two vertical brown lines indicate 
when the numerical simulation reaches mui = 0.063 and 0.1, 
respectively; the labels along the top horizontal axes give the 
number of gravitational- wave cycles before muj = 0.1. 
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B. Comparing different post-Newtonian 
approximants 

The previous section presented detailed comparisons 
of our numerical waveforms with four different post- 
Newtonian approximants. We now turn our attention to 
some comparisons between these approximants. In this 
section we also explore further how the post-Newtonian 
order influences agreement between numerical and post- 
Newtonian waveforms. 

Figure [221 presents phase differences as a function of 
time for all four PN approximants we consider here and 
for different PN orders. The post-Newtonian and nu- 
merical waveforms are matched at mujm = 0.04, about 
9 cycles after the beginning of the numerical waveform, 
and about 21 cycles before its end. We find that some 
PN approximants at some particular orders agree exceed- 
ingly well with the numerical results. The best match 
is easily TaylorT4 at 3.5PN order, and the next best 
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To get the complete waveform to 3PN order, only the (2, 2) mode 
must be known to 3PN order; other modes must be known to 
smaller PN orders. For an equal mass, non-spinning binary, all 
modes except the (3, 2) mode are currently known to sufficient 
order to get a complete 3PN waveform |62|| . 



FIG. 23: Same as Fig. 1221 but showing only the last stage of 
the inspiral. The horizontal axis ends at the estimated time 
of merger, {t — r*)cAH = 3950m, cf. Sec 111 Gl The top and 
bottom panels use different vertical scales. 
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match is Taylor T4 at 2.0PN order. Some approximants 
behave significantly worse, such as the TaylorTl and Tay- 
lorT4 waveforms at 2.5PN order. The 2.5PN and 3PN 
TaylorT3 waveforms agree very well with the numerical 
waveform at early times, but at late times they accumu- 
late a large phase difference; the 2.5PN TaylorTS wave- 
form ends even before the numerical waveform reaches 
muj = 0.1 (the rightmost vertical brown line in Fig. I22p. 

We also find that all four PN approximants, when com- 
puted to 3PN order or higher, match the numerical wave- 
form (and each other) quite closely at early times, when 
all PN approximants are expected to be accurate. How- 
ever, at late times, t — r*> 2500to, the four PN approx- 
imants begin to diverge, indicating that PN is beginning 
to break down. 

Figure [23] is an enlargement of Fig. [22] for the last 
10 gravitational wave cycles before merger. This figure 
shows in more detail how the different PN approximants 
behave near merger. 

Figure [Ml presents similar results in a different format. 
We compute the phase differences between the numerical 
waveform and the various post-Newtonian approximants 
at the times when the numerical waveform reaches grav- 
itational wave frequencies muj — 0.063 and muj — 0.1 
(the times corresponding to these frequencies are also in- 
dicated by brown lines in Fig. [52]). We then plot these 
phase differences as a function of the post-Newtonian or- 
der (using equal order in phase and amplitude, except for 
3.5PN order, where we use 3.0PN in amplitude). Three 
PN approximants end before to.i- TaylorTl 2.0/2.0, Tay- 
lorT3 2.5/2.5, TaylorT3 3.5/3.0. These data points there- 
fore cannot be included in the right panel of Fig. [24] 

The general trend seen in Fig. [24] is that the phase dif- 
ference decreases with increasing PN order. However, 
this convergence is not monotonic, and the scatter in 
Fig. [24] can be larger than the phase differences them- 
selves. For example, the OPN waveforms are about as 
good as the 2.5PN waveforms for TaylorTl and TaylorT4, 
and the 2PN TaylorT4 waveform agrees with the numer- 
ical results much better than do either the 2.5PN or 3PN 
TaylorT4 waveforms. Considering Fig. [M] it seems diffi- 
cult to make statements about the convergence with PN 
order for any particular PN approximant, or statements 
about which PN orders are generally "good" . Given that 
at fixed PN order the different approximants differ merely 
by the treatment of uncontrolled higher-order terms, the 
scatter in Fig. [M] in some sense represents the truncation 
error at each PN order. While some PN approximants 
at certain orders may show better agreement with the 
numerical simulation, we are not aware of any means to 
predict this besides direct comparisons to numerical sim- 
ulations (as is done here). In particular. Fig. [24l suggests 
that the remarkable agreement between our numerical re- 
sults and the 3.5PN TaylorT4 approximant may be sim- 
ply due to luck; clearly, more PN-NR comparisons are 
needed, with different mass ratios and spins, to see if 
this is the case. 
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FIG. 24: Phase differences between numerical and post- 
Newtonian waveforms at two selected times close to merger. 
Waveforms are matched at muJm = 0.04, and phase differ- 
ences are computed at the time when the numerical simulation 
reaches muj = 0.063 (left panel) and muj = 0.1 (right panel). 
Differences are plotted versus PN order (equal order in phase 
and amplitude, except the '3.5 PN' points are 3.5/3.0). On 
the right plot, the IPN data points are off scale, clustering at 
— 15 radians. The thin black bands indicate upper bounds on 
the uncertainty of the comparison as discussed in Sec. IV Al 



VII. CONCLUSIONS 

We have described numerical simulations of an equal 
mass, non-spinning binary black hole spacetime cover- 
ing 15 orbits of inspiral just prior to the merger of the 
two black holes. Using a multi-domain pseudospectral 
method we are able to extract the gravitational wave 
content measured by a distant observer with a phase ac- 
curacy of better than 0.02 radians over the roughly 30 cy- 
cles of gravitational radiation observed. We demonstrate 
that in order to achieve this accuracy it is necessary to 
accurately extrapolate the waveform from data obtained 
at extraction surfaces sufficiently far from the center of 
mass of the system. When comparing to zero-spin, zero- 
eccentricity PN formulae, our phase uncertainty increases 
to 0.05 radians because the numerical simulation has 
a small but nonzero orbital eccentricity and small but 
nonzero spins on the holes. 

Judging from the case in which we match at mujm = 
0.04, our numerical simulations are consistent (within our 
estimated phase uncertainty) with all PN approximants 
(at the highest PN order) from the beginning of our inspi- 
ral until about 15 gravitational wave cycles prior to the 
merger of the binary. This agreement provides an impor- 
tant validation of our numerical simulation. It also estab- 
lishes a regime in which the 3.5-th order post-Newtonian 
waveforms are accurate to this level, at least for an equal 
mass, non-spinning black hole binary. After this point, 
the various PN approximants begin to diverge, suggest- 
ing that the approximation is beginning to break down. 
Since there are many different PN approximants (includ- 
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ing Fade ^ and effective-one-body [H, [H [H, ill which 
were not discussed in this paper) it may be possible to 
find a clever way to push the PN expansion beyond its 
breaking point. 

Indeed, we find that one approximant, TaylorT4 at 
3.5PN in phase, works astonishingly well, agreeing with 
our numerical waveforms for almost the entire 30-cycle 
length of our runs. Given the wide scatter plot of pre- 
dictions by various PN approximants, it is likely that 
TaylorT4 3.5/3.0 simply got lucky for the equal mass 
non-spinning black hole binary. In fact, the assumption 
of adiabaticity (i.e., circular orbits) is known to lead to 
much larger phase differences r elati ve to a non-adiabatic 
inspiral (see Fig. 4 of [HI] and [l43l |) than the phase dif- 
ferences between NR and TaylorT4 wc find in Fig. \T9\ 
Thus it seems that the uncontrolled higher order terms 
of TaylorT4 3.5/3.0 balance the error introduced by the 
adiabaticity assumption to a remarkable degree. It re- 
mains to be determined whether similar cancellations oc- 
cur when the black hole masses are unequal or when the 
holes have nonzero spin. 

Regardless of the robustness of TaylorT4, it seems ev- 
ident that numerical simulations are needed in order to 
know which, if any, PN approximant yields the correct 
waveform after the various approximants begin to di- 
verge. For there is no a priori reason why Taylor T4 
should be a better choice than TaylorTl as they differ 
only in whether the ratio of gravitational wave flux to 
the derivative of the orbital energy with respect to fre- 
quency is left as a ratio of post-Newtonian expansions or 
re-expanded as a single post-Newtonian expansion. 

The surprising accuracy of TaylorT4 3.5/3.0 in the 
gravitational frequency range from muj — 0.035 through 
mu) = 0.15, for the equal mass, non-spinning inspiral of 
two black holes, in principle could form a basis for eval- 
uating the errors of numerical simulations. Instead of 
worrying about errors due to different formulations, ini- 
tial data, boundary conditions, extraction methods, etc., 
perhaps a long inspiral simulation could be compared 
with TaylorT4 3.5/3.0 in order to get a direct estimate 
of the phase error. Similarly, because of its good agree- 
ment, TaylorT4 3.5/3.0 could also be used to address 
questions that require much longer waveforms than cur- 
rently available, for instance the question of when lower 
order post-Newtonian waveforms become unreliable. 

We find that the 3PN contributions to the amplitude 
of the (2, 2) modes improve their accuracy with respect to 
the numerical waveforms. This suggests that for accurate 
parameter estimation, it may be desirable to compute 
the full 3PN amplitude for the polarization waveforms. 
Despite the formidable nature of the calculation required, 
it would also be interesting to see how the inclusion of 
4PN order corrections to the phasing would affect our 
comparisons. 

Much work still needs to be done to improve the com- 
parison between NR and PN. Our primary goal is to push 
our simulations through merger and ringdown so that we 



may compare various resummed PN approximants and 
the effective-one-body approximants during the last cy- 
cle of inspiral and merger, as well as test Taylor T4 3.5/3.0 
closer to merger. We also intend to do long inspirals with 
arbitrary masses and spins in order to test the robustness 
of PN over a range of these parameters. 

Furthermore we wish to improve our initial data. 
There is a large amount of 'junk radiation' present in the 
initial data that limits how early we can match PN and 
NR waveforms. Reduction of this junk radiation jl44l | 
would improve the accuracy of our simulations as well. 

Finally, we have done just a simple comparison be- 
tween NR and PN, without including any treatment of 
effects that are important for real gravitational wave de- 
tectors such as limited bandwidth and detector noise. In 
order to more directly address the suitability of PN for- 
mulae for analyzing data from gravitational wave detec- 
tors, it will be necessary to fold in the properties of the 
detector, to consider specific values for the total mass of 
the binary, and to fit for the mass from the waveforms 
rather than assuming that the PN and NR waveforms 
correspond to the same mass. We leave this for future 
work. 
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